Attachment 'getlandaupeak.C'

Download

void getlandaupeak(TH2* hist, int idsta=1, TString opt=""){
  /*
    ROOT Macro for obtaining Landau peak values from 2d histogram (id vs Qav)
    input: 2d-hist
           option="nstop" -> wait is skipped
    output: landaupeak.txt  -> peak values
  */

  Int_t nbin=hist->GetXaxis()->GetNbins();
  Int_t num=0;
  vector<double> peak;//Landau peak value
  if (nbin>1000){
    cout<<"Too much X-axis bins. Please check the input 2D-histogram."<<endl;
    exit(0);
  }

  for (Int_t i=idsta-1;i<nbin;i++){// begin of loop ------->
    num++;
    Int_t id = i+1;
    TH1 *h=MakeBanY(hist,id);

    Double_t p[3];
    FitHist(h,p);
    peak.push_back(p[1]);

    // Draw
    h->Draw();
    gPad->Update();
    cout<<"ID="<<id<<" peak="<<p[1]<<endl;

    if (opt!="nstop")
      if(Wait()) break;

  }// end of loop <--------------

  // file write
  ofstream ofs1;
  ofs1.open("landaupeak.txt");
  vector<double>::iterator it=peak.begin();
  while(it!=peak.end()){
    ofs1<<*it<<endl;
    it++;
  }
  ofs1.close();

  cout<<num<<" peaks are fitted, END"<<endl;
}
//___________________________________________________________________________
void FitHist(TH1 *h, double p[]){

  int searchedbin[2];
  SearchBins(h,searchedbin);// search local max. and min.

  TF1 *func = new TF1("myfunc","landau");
  double mean = h->GetMean();
  double rms = h->GetRMS();
  func->SetParameters( h->GetMaximum(),
                       h->GetBinCenter(searchedbin[0]),
                       h->GetRMS()      );
  // fitting range
  Double_t xmin,xmax;
  xmin = h->GetBinCenter(searchedbin[1]);
  xmax = h->GetBinCenter(searchedbin[0])+rms;

  h->Fit(func,"Q","",xmin,xmax);
  func->GetParameters(p);
}
//___________________________________________________________________________
void SearchBins(TH1* h, int searchedbin[]){
  /*
    search bin at local maximum (Landau peak) and local minimum 
    from high energy side
  */
  const int range=20;  // for taking average of contents
  int thre = h->GetMaximum()*0.5;// Threshold value 

  Int_t nbin=h->GetXaxis()->GetNbins();
  int isearch=0;
  int sum,sumold=0;
  for(int i=nbin;i>range;i--){
    if (h->GetBinContent(i)<thre) continue;
    sum=h->Integral(i-range,i);
    if (isearch==0 && sum<sumold){// Local maximum
      searchedbin[0]=i;
      isearch=1;
    } else if (isearch==1 && sum>sumold){// Local minimum
      searchedbin[1]=i;
      break;
    }
    sumold=sum;
  }
  if (isearch==0){
    searchedbin[0]=h->GetMaximumBin();
    searchedbin[1]=0;
  }
}
//___________________________________________________________________________
TH1* MakeBanY(TH2* hist, Int_t id){
  TString str=hist->GetName();
  ostringstream oss;
  oss<<id;
  str+="_bny"+oss.str();
  return hist->ProjectionY(str.Data(),id,id);
}
//___________________________________________________________________________
bool Wait(){
  //check and quit if q is typed
  printf("[enter]:continue, [q]:quit >");
  char s[1];
  gets(s);
  if(s[0]=='q') {
    cout<<"Break fit loop"<<endl;
    return true;
  }
  return false;
}
//___________________________________________________________________________

Attached Files

To refer to attachments on a page, use attachment:filename, as shown below in the list of files. Do NOT use the URL of the [get] link, since this is subject to change and can break easily.

• [get | view] (2012-11-02 05:21:15, 1.6 KB) [[attachment:BPCHist.C]]
• [get | view] (2013-02-27 09:18:59, 0.8 KB) [[attachment:CoinHist.C]]
• [get | view] (2012-11-02 05:10:49, 2.6 KB) [[attachment:FDC2Hist.C]]
• [get | view] (2012-11-02 04:31:49, 2.1 KB) [[attachment:HODHist.C]]
• [get | view] (2012-11-02 04:31:32, 1.6 KB) [[attachment:ICBHist.C]]
• [get | view] (2014-12-31 11:41:13, 19.5 KB) [[attachment:OnlineMonitor.cc]]
• [get | view] (2014-12-31 11:41:02, 2.6 KB) [[attachment:OnlineMonitor.hh]]
• [get | view] (2012-11-02 05:25:41, 1.7 KB) [[attachment:PlaHist.C]]
• [get | view] (2012-02-15 09:14:24, 3.2 KB) [[attachment:ch2ns.C]]
• [get | view] (2012-02-15 10:05:13, 3.0 KB) [[attachment:getlandaupeak.C]]
• [get | view] (2012-02-15 09:14:53, 2.7 KB) [[attachment:getped.C]]
• [get | view] (2012-02-13 09:49:04, 2.6 KB) [[attachment:makevsta.C]]
• [get | view] (2012-02-19 14:44:04, 0.5 KB) [[attachment:slew.C]]

You are not allowed to attach a file to this page.