EICMCAnalysis/analysis_eicsmear/eicsmear_check_calorimeter.C

0001 int
0002 eicsmear_check_calorimeter( TString filename_mc,
0003                 TString filename_mc_smeared = "")
0004 {
0005   /* Uncomment this line when running without compilation */
0006   gSystem->Load("libeicsmear");
0007
0008   /* Open input files. */
0009   TFile *file_mc = new TFile(filename_mc, "OPEN");
0010   TFile *file_mc_smeared = new TFile(filename_mc_smeared, "OPEN");
0011
0012   /* Get trees from files. */
0013   TTree *tree = (TTree*)file_mc->Get("EICTree");
0014   TTree *tree_smeared = (TTree*)file_mc_smeared->Get("Smeared");
0015
0016   /* Add friend to match branches in trees. */
0017   tree->AddFriend(tree_smeared);
0018
0019   erhic::EventPythia *event  = NULL;
0020   Smear::Event       *eventS = NULL;
0021
0022   tree->SetBranchAddress("event", &event);
0023   tree->SetBranchAddress("eventS", &eventS);
0024
0025   /* Define base cut for reconstructed final state particles */
0026   TCut cut_base("particles.KS==1 && Smeared.particles.p > 0 && (particles.id==11 || particles.id==22)");
0027
0028   /* colors array */
0029   unsigned colors[7] = {1,2,3,4,6,7,14};
0030
0031   /* select how many eta settings to plot */
0032   unsigned etas_plotmax = 7;
0033
0034   /* Create vector of theta values to include */
0035   vector< double > etas;
0036   etas.push_back(-2.75);
0037   etas.push_back(-2.25);
0038   etas.push_back(-1.75);
0039   etas.push_back(-0.25);
0040   etas.push_back( 0.25);
0041   etas.push_back( 1.75);
0042   etas.push_back( 2.25);
0043
0044   /* Create fit function */
0045   TF1* f_momres = new TF1("f_momres", "sqrt( ([0])**2 + ([1]/(x**0.5))**2 )" );
0046
0047   cout << "\nFit function: " << f_momres->GetTitle() << "\n" << endl;
0048
0049   /* Create scratch canvas */
0050   TCanvas *cscratch = new TCanvas("cscratch");
0051
0052   /* Create framehistogram */
0053   TH1F* hframe = new TH1F("hframe","",100,0,40);
0054   hframe->GetYaxis()->SetRangeUser(0,0.15);
0055   hframe->GetXaxis()->SetTitle("Energy (GeV/c)");
0056   hframe->GetYaxis()->SetTitle("#sigma_{E}/E");
0057
0058   /* Create 2Dhistogram */
0059   unsigned nbinsp = 40;
0060   TH2F* h2tmp = new TH2F("h2tmp","",nbinsp,0,(float)nbinsp,50,-1,1);
0061   h2tmp->GetXaxis()->SetTitle("Energy (GeV)");
0062   h2tmp->GetYaxis()->SetTitle("(#Delta E)/E");
0063
0064   /* create combined canvas plot */
0065   TCanvas *c1 = new TCanvas();
0066   hframe->Draw();
0067
0068   /* Create legend */
0069   TLegend* leg_eta = new TLegend( 0.25, 0.40, 0.45, 0.90);
0070
0071   /* Create resolution-vs-momentum plot with fits for each selected theta value */
0072   for ( int i = 0; i < etas.size(); i++ )
0073     {
0074       /* Switch to scratch canvas */
0075       cscratch->cd();
0076
0077       double eta = etas.at(i);
0078
0079       cout << "\n***Eta = " << eta << endl;
0080
0081       /* Define range of theta because float comparison with fixed value doesn't work
0082          too well for cuts in ROOT trees */
0083       double eta_min = eta - 0.25;
0084       double eta_max = eta + 0.25;
0085
0086       /* Cut for tree */
0087       TCut cutx( Form("particles.p > 1 && ( particles.eta > %f && particles.eta < %f )", eta_min, eta_max) );
0088       cout << cutx.GetTitle() << endl;
0089
0090       /* Fill 2D histogram and fit slices with gaussian to get track resolution */
0091       h2tmp->Reset();
0092       tree->Draw("(Smeared.particles.E-particles.E)/particles.E:particles.E >> h2tmp",cut_base && cutx,"");
0093       h2tmp->FitSlicesY();
0094
0095       /* Get histogram with sigma from gauss fits */
0096       TH1D* h2tmp_2 = (TH1D*)gDirectory->Get("h2tmp_2");
0097
0098       /* Fill fit parameters in TGraphErrors */
0099       TGraphErrors *gres = new TGraphErrors(nbinsp);
0100       gres->SetMarkerColor(colors[i]);
0101       for ( unsigned bini = 1; bini < nbinsp; bini++ )
0102     {
0103       double sigm_i =     h2tmp_2->GetBinContent(bini);
0104       double sigm_err_i = h2tmp_2->GetBinError(bini);
0105
0106       gres->SetPoint( bini-1, h2tmp->GetXaxis()->GetBinCenter(bini), sigm_i );
0107       gres->SetPointError( bini-1, 0, sigm_err_i );
0108     }
0109
0110       /* Only plot first few lines; if not plotting, still do the fit */
0111       if ( i < etas_plotmax )
0112         {
0113           /* Add graph to legend */
0114           leg_eta->AddEntry(gres, Form("#eta = %.1f", eta), "P");
0115
0116           /* Add graph to plot */
0117           c1->cd();
0118           gres->Draw("Psame");
0119           f_momres->SetLineColor(colors[i]);
0120           gres->Fit(f_momres);
0121         }
0122       else
0123         {
0124           gres->Fit(f_momres);
0125         }
0126     }
0127
0128   /* Draw legend */
0129   TCanvas *c2 = new TCanvas();
0130   leg_eta->Draw();
0131
0132   return 0;
0133 }