analysis/FrawleyTrackEvalMacros/project_quarkonium_mass.C

0001 void project_quarkonium_mass()
0002 {
0003   //SetsPhenixStyle();
0004   gROOT->SetStyle("Plain");
0005   gStyle->SetOptStat(0);
0006   gStyle->SetOptFit(0);
0007   gStyle->SetOptTitle(0);
0008
0009   TFile *fout = new TFile("root_files/ntp_quarkonium_out.root","recreate");
0010
0011   TH1D* recomass[4];
0012   char name[500];
0013   for(int istate = 0; istate < 4; ++istate)
0014     {
0015       sprintf(name, "recomass%i",istate);
0016       //recomass[istate] = new TH1D(name,name,100,7.0,11.0);
0017       recomass[istate] = new TH1D(name,name,200,7.0,11.0);
0018
0019     }
0020
0021   TH2D *hdca2d = new TH2D("hdca2d","hdca2d", 200, 0.0, 50.0, 1000, -0.002, 0.002);
0022   TH2D *hpcaz = new TH2D("hpcaz","hpcaz",200, 0.0, 50.0, 1000,-0.1,0.1);
0023
0024   Float_t event;
0025   Float_t px;
0026   Float_t py;
0027   Float_t pz;
0028   Float_t pt;
0029   Float_t gembed;
0030   Float_t quality;
0031   Float_t charge;
0032   Float_t dca2d;
0033   Float_t pcaz;
0034   Float_t gvz;
0035   Float_t ntpc;
0036   Float_t nmaps;
0037
0038   Float_t decaymass=0.000511;
0039
0040   // what do we have in the file?
0041   int events_per_file = 1;
0042   int num_ups_per_event = 20;
0043
0044  // what maximum yields do we want?
0045   // These are pp yields for 197 /nb of pp running, with
0046   //double pair_acc = 0.60 * 0.9 * 0.98;  // fraction into 1 unit of rapidity * pair eID fraction * trigger fraction
0047
0048   // pp case:
0049   // These are pp yields from Sasha's TN for 13,100 B pp collisions
0050   //int Npp[3] = {17074, 4440, 2475};  // 1S, 2S, 3S
0051
0052   // AuAu case:
0053   // These are yields from Sasha's TN for 24B AuAu events with 100% eID
0054   int Npp[3] = {27120, 7053, 3932};  // 1S, 2S, 3S unsuppressed
0055
0056   double suppr_state[3] = {0.535, 0.17, 0.035}; // suppression factor from Strickland paper
0057   double eID_eff = 0.9;
0058   double tracking_eff = 0.95;
0059
0060   double quality_cut = 10;
0061
0062   bool no_suppression = true;
0063
0064   int nupsmax[3] = {0,0,0};
0065   for(int istate = 0; istate < 3; ++istate)
0066     {
0067       if(no_suppression) suppr_state[istate] = 1.0;
0068
0069       nupsmax[istate] = eID_eff * tracking_eff * suppr_state[istate] * Npp[istate];
0070
0071       nupsmax[istate] = 100000;  // take all
0072       cout << "Requested max Upsilon counts for istate " << istate << " = " << nupsmax[istate] << endl;
0073     }
0074
0075
0076   int nups[3] = {0,0,0};
0077
0078   cout << "Upsilons requested = " << nupsmax << endl;
0079
0080   for(int ifile=0;ifile<1000;ifile++)
0081     {
0082       if(nups[0] > nupsmax[0])
0083     break;
0084
0085       char name[2000];
0086
0087       //sprintf(name,"/sphenix/user/frawley/cluster_efficiency/macros/macros/g4simulations/eval_output/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
0088       //sprintf(name,"/sphenix/user/frawley/acts_qa/macros/macros/g4simulations/eval_output_3/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
0089       sprintf(name,"/sphenix/user/frawley/acts_qa/macros/macros/g4simulations/eval_output_2/g4svtx_eval_%i.root_g4svtx_eval.root",ifile);
0090
0091       cout << "Adding file number " << ifile << " with name " << name << endl;
0092
0093       TChain* ntp_vertex = new TChain("ntp_vertex","reco events");
0094       ntp_vertex->Add(name);
0095
0096       TChain* ntp_track = new TChain("ntp_track","reco tracks");
0097       ntp_track->Add(name);
0098
0099       ntp_track->SetBranchAddress("event",&event);
0100       ntp_track->SetBranchAddress("px",&px);
0101       ntp_track->SetBranchAddress("py",&py);
0102       ntp_track->SetBranchAddress("pz",&pz);
0103       ntp_track->SetBranchAddress("pt",&pt);
0104       ntp_track->SetBranchAddress("charge",&charge);
0105       ntp_track->SetBranchAddress("dca2d",&dca2d);
0106       ntp_track->SetBranchAddress("pcaz",&pcaz);
0107       ntp_track->SetBranchAddress("gvz",&gvz);
0108       ntp_track->SetBranchAddress("gembed",&gembed);
0109       ntp_track->SetBranchAddress("quality",&quality);
0110       ntp_track->SetBranchAddress("ntpc",&ntpc);
0111       ntp_track->SetBranchAddress("nmaps",&nmaps);
0112
0113       int ntracks=ntp_track->GetEntries();
0114
0115       // capture the dca2d and dcaz histograms for the pions
0116       for(int i=0;i<ntracks;i++)
0117     {
0118       ntp_track->GetEntry(i);
0119
0120       if(gembed != 2)
0121         continue;
0122
0123       hdca2d->Fill(pt, dca2d);
0124       hpcaz->Fill(pt, pcaz-gvz);
0125     }
0126
0127       //events_per_file = ntp_vertex->GetEntries();
0128       cout << " events in this file = " << events_per_file << endl;
0129
0130       for(int iev = 0;iev<events_per_file;iev++)
0131     {
0132       for(int iups = 0;iups< num_ups_per_event; iups++)
0133         {
0134           //int istate = iups%3;
0135           int istate = 0;
0136
0137           if(nups[istate] > nupsmax[istate])
0138         continue;
0139
0140           int nlept1 = 0;
0141           int nlept2 = 0;
0142           int lept1=-1;
0143           int lept2 = -1;
0144           Float_t px1=0;
0145           Float_t px2=0;
0146           Float_t py1=0;
0147           Float_t py2=0;
0148           Float_t pz1 = 0;
0149           Float_t pz2 = 0;
0150
0151           TLorentzVector t1;
0152           TLorentzVector t2;
0153           for(int i=0;i<ntracks;i++)
0154         {
0155           ntp_track->GetEntry(i);
0156
0157           if(event != iev)
0158             continue;
0159
0160           if(gembed != iups+3)
0161             continue;
0162
0163           if(quality > quality_cut)
0164             continue;
0165
0166           if(ntpc < 20)
0167             continue;
0168
0169           if(nmaps <= 1)
0170             continue;
0171
0172           if(sqrt(px*px + py*py + pz*pz) < 1.0)
0173             continue;
0174
0175           //           cout << "event = " << event << " iups " << iups << " charge " << charge << " px " << px << " py " << py << " pz " << pz << " gembed " << gembed << " quality " << quality << endl;
0176
0177
0178           if(charge == 1)
0179             {
0180               nlept1++;
0181
0182               lept1 = 1;
0183               px1 = px;
0184               py1 = py;
0185               pz1 = pz;
0186
0187               Float_t E1 = sqrt( pow(px1,2) + pow(py1,2) + pow(pz1,2) + pow(decaymass,2));
0188               t1.SetPxPyPzE(px1,py1,pz1,E1);
0189             }
0190
0191           if(charge == -1)
0192             {
0193               nlept2++;
0194
0195               lept2 = 1;
0196               px2 = px;
0197               py2 = py;
0198               pz2 = pz;
0199
0200               Float_t E2 = sqrt( pow(px2,2) + pow(py2,2) + pow(pz2,2) + pow(decaymass,2));
0201               t2.SetPxPyPzE(px2,py2,pz2,E2);
0202             }
0203           //cout << " lept1 " << lept1 << " lept2 " << lept2 << endl;
0204         }
0205
0206           // calculate mass
0207           if(nlept1 == 1 && nlept2 == 1)
0208           //if(nlept1 > 0 && nlept2 > 0)
0209         {
0210           TLorentzVector tsum;
0211           tsum = t1+t2;
0212           recomass[istate]->Fill(tsum.M());
0213           recomass[3]->Fill(tsum.M()); // total spectrum
0214           if(tsum.M() > 7 && tsum.M() < 11)
0215             nups[istate]++;
0216           //cout << " event " << event << " gembed " << gembed << endl;
0217           //cout << " istate " << " nups[istate] " << nups[istate] << " reco mass " << tsum.M() << endl;
0218         }
0219         }
0220     }
0221       delete ntp_track;
0222       delete ntp_vertex;
0223     }
0224
0225   cout << " requested: " << nupsmax[0] << " reconstructed: " << nups[0] << " Upsilons " << endl;
0226
0227   for(int istate =0; istate < 3; ++istate)
0228     {
0229       recomass[istate]->Write();
0230       cout << " requested for istate: " << istate << " nupsmax = " << nupsmax[istate] << " reconstructed = " << nups[istate] << " Upsilons " << endl;
0231     }
0232   recomass[3]->Write();
0233
0234   hdca2d->Write();
0235   hpcaz->Write();
0236
0237   fout->Close();
0238
0239 }
0240
0241