macros/tracking_analysis/TrackParametrization.C

0001 int
0002 TrackParametrization( TString csvfile="fitslices_out.csv" )
0003 {
0004
0005   /* Read data from input file */
0006   TTree *tres = new TTree();
0007   tres->ReadFile( csvfile, "ptrue:etatrue:psig:psig_err:pmean:pmean_err:norm", ',' );
0008
0009   /* Print read-in tree */
0010   tres->Print();
0011
0012   /* colors array */
0013   unsigned colors[8] = {1,2,3,4,6,7,14,16};
0014
0015   /* Create vector of theta values to include for visualization*/
0016   vector< double > etas_vis;
0017   etas_vis.push_back(-2.75);
0018   etas_vis.push_back(-2.25);
0019   etas_vis.push_back(-1.75);
0020   etas_vis.push_back(-0.25);
0021   etas_vis.push_back( 0.25);
0022   etas_vis.push_back( 1.75);
0023   etas_vis.push_back( 2.25);
0024
0025 //  etas_vis.push_back(-3.25);
0026 //  etas_vis.push_back(-2.25);
0027 //  etas_vis.push_back(-1.25);
0028 //  etas_vis.push_back(-0.25);
0029 //  etas_vis.push_back( 0.25);
0030 //  etas_vis.push_back( 1.25);
0031 //  etas_vis.push_back( 2.25);
0032 //  etas_vis.push_back( 3.25);
0033
0034   /* Create vector of theta values to include for fitting*/
0035   vector< double > etas_fit;
0036   for ( double eta = -4.45; eta < 4.5; eta += 0.1 )
0037     etas_fit.push_back( eta );
0038
0039   /* Create fit function */
0040   TF1* f_momres = new TF1("f_momres", "sqrt( [0]*[0] + [1]*[1]*x*x )" );
0041
0042   cout << "\nFit function: " << f_momres->GetTitle() << "\n" << endl;
0043
0044   /* Create scratch canvas */
0045   TCanvas *cscratch = new TCanvas("cscratch");
0046
0047   /* Create framehistogram */
0048   TH1F* hframe = new TH1F("hframe","",100,0,40);
0049   hframe->GetYaxis()->SetRangeUser(0,0.15);
0050   hframe->GetYaxis()->SetNdivisions(505);
0051   hframe->GetXaxis()->SetTitle("Momentum (GeV/c)");
0052   hframe->GetYaxis()->SetTitle("#sigma_{p}/p");
0053
0054   /* create combined canvas plot */
0055   TCanvas *c1 = new TCanvas();
0056   hframe->Draw();
0057
0058   /* Create legend */
0059   TLegend* leg_eta = new TLegend( 0.2, 0.6, 0.5, 0.9);
0060   leg_eta->SetNColumns(2);
0061
0062   /* Create ofstream to write fit parameter results */
0063   ofstream ofsfit("track_momres_new.csv");
0064   ofsfit<<"eta,par1,par1err,par2,par2err"<<endl;
0065
0066   /* Create resolution-vs-momentum plot with fits for each selected theta value */
0067   for ( int i = 0; i < etas_fit.size(); i++ )
0068     {
0069       /* Switch to scratch canvas */
0070       cscratch->cd();
0071
0072       double eta = etas_fit.at(i);
0073
0074       /* No tracking outside -4 < eta < 4 */
0075       if ( eta < -4 || eta > 4 )
0076     continue;
0077
0078       cout << "\n***Eta = " << eta << endl;
0079
0080       /* Define range of theta because float comparison with fixed value doesn't work
0081      too well for cuts in ROOT trees */
0082       double eta_min = eta * 0.999;
0083       double eta_max = eta * 1.001;
0084
0085       /* Cut for tree */
0086       TCut cutx( Form("ptrue > 1 && ( (etatrue > 0 && (etatrue > %f && etatrue < %f)) || (etatrue < 0 && (etatrue < %f && etatrue > %f)) )", eta_min, eta_max, eta_min, eta_max) );
0087
0088       /* "Draw" tree on scratch canvas to fill V1...V4 arrays */
0089       tres->Draw("psig:ptrue:psig_err:0", cutx );
0090
0091       /* Create TGraphErrors with selected data from tree */
0092       TGraphErrors *gres = new TGraphErrors( tres->GetEntries(cutx),
0093                          &(tres->GetV2())[0],
0094                          &(tres->GetV1())[0],
0095                          &(tres->GetV4())[0],
0096                          &(tres->GetV3())[0] );
0097
0098       /* reset function parameters before fit */
0099       f_momres->SetParameter(0,0.1);
0100       f_momres->SetParameter(1,0.1);
0101
0102       /* Only plot pseudorapidities listed on etas_vis; if not plotting, still do the fit */
0103       bool vis = false;
0104       int vi = 0;
0105
0106       for ( vi = 0; vi < etas_vis.size(); vi++ )
0107     {
0108       if ( abs( etas_vis.at(vi) - eta ) < 0.001 )
0109         {
0110           vis = true;
0111           break;
0112         }
0113     }
0114
0115       if ( vis )
0116     {
0117       /* Add graph to legend */
0118       leg_eta->AddEntry(gres, Form("#eta = %.1f", eta), "P");
0119
0120       /* Add graph to plot */
0121       c1->cd();
0122       gres->SetMarkerColor(colors[vi]);
0123       gres->Draw("Psame");
0124       f_momres->SetLineColor(colors[vi]);
0125       gres->Fit(f_momres);
0126     }
0127       else
0128     {
0129       gres->Fit(f_momres);
0130     }
0131
0132       /* Write fir results to file */
0133       double par1 = f_momres->GetParameter(0);
0134       double par1err = f_momres->GetParError(0);
0135       double par2 = f_momres->GetParameter(1);
0136       double par2err = f_momres->GetParError(1);
0137       ofsfit << eta << "," << par1 << "," << par1err << "," << par2 << "," << par2err << endl;
0138
0139     }
0140
0141   /* Draw legend */
0142   c1->cd();
0143   //TCanvas *c2 = new TCanvas();
0144   //hframe->Draw();
0145   leg_eta->Draw();
0146
0147   /* Print plots */
0148   c1->Print("track_momres_vareta.eps");
0149   //c2->Print("track_momres_vareta_legend.eps");
0150
0151   /* Close output stream */
0152   ofsfit.close();
0153
0154   return 0;
0155 }