File indexing completed on 2025-08-05 08:12:06
0001 #include <vector>
0002 int eicsmear_dvmp_tree(TString filename_mc = "/sphenix/user/gregtom3/data/Fall2018/JPsi_reco_studies/reconstructedQ2/18x275_DVMP_1M_ascii_converted.root",
0003 TString filename_mc_smeared = "/sphenix/user/gregtom3/data/Fall2018/JPsi_reco_studies/reconstructedQ2/18x275_DVMP_1M_ascii_converted.smear.root")
0004 {
0005
0006
0007
0008
0009 gSystem->Load("libeicsmear");
0010 gROOT->LoadMacro("./sPHENIXStyle/sPhenixStyle.C");
0011 SetsPhenixStyle();
0012
0013
0014
0015
0016
0017
0018
0019
0020 TFile *file_mc = new TFile(filename_mc, "OPEN");
0021 TFile *file_mc_smeared = new TFile(filename_mc_smeared, "OPEN");
0022
0023
0024 TTree *tree = (TTree*)file_mc->Get("EICTree");
0025 Double32_t Q2fromEICTree;
0026 tree->SetBranchAddress("trueQ2",&Q2fromEICTree);
0027 TTree *tree_smeared= (TTree*)file_mc_smeared->Get("Smeared");
0028
0029
0030 TFile *myfile = new TFile("/sphenix/user/gregtom3/data/Fall2018/JPsi_reco_studies/reconstructedQ2/analysisTree.root","RECREATE");
0031
0032 tree->AddFriend(tree_smeared);
0033
0034
0035
0036 TTree *theTree = new TTree("DVMP_Tree","");
0037 Float_t de_phi_truth, de_eta_truth, de_pt_truth, de_p_truth;
0038 Float_t dp_phi_truth, dp_eta_truth, dp_pt_truth, dp_p_truth;
0039 Float_t sp_phi_truth, sp_eta_truth, sp_pt_truth, sp_p_truth;
0040 Float_t se_phi_truth, se_eta_truth, se_pt_truth, se_p_truth;
0041 Float_t jpsi_phi_truth, jpsi_eta_truth, jpsi_pt_truth, jpsi_p_truth;
0042 Float_t de_phi_reco, de_eta_reco, de_pt_reco, de_p_reco;
0043 Float_t dp_phi_reco, dp_eta_reco, dp_pt_reco, dp_p_reco;
0044 Float_t sp_phi_reco, sp_eta_reco, sp_pt_reco, sp_p_reco;
0045 Float_t se_phi_reco, se_eta_reco, se_pt_reco, se_p_reco;
0046 Float_t jpsi_phi_reco, jpsi_eta_reco, jpsi_pt_reco, jpsi_p_reco;
0047 Float_t invariant_reco_decay;
0048 Float_t invariant_reco_scatter;
0049 Float_t jpsi_px_truth, jpsi_py_truth, jpsi_pz_truth;
0050 Float_t jpsi_px_reco, jpsi_py_reco, jpsi_pz_reco;
0051 Double_t Q2;
0052
0053 Float_t t;
0054 theTree->Branch("Q2",&Q2,"Event QSquared");
0055 theTree->Branch("t",&t,"Event t");
0056 theTree->Branch("de_phi_truth",&de_phi_truth,"Decay Electron Truth Phi");
0057 theTree->Branch("de_eta_truth",&de_eta_truth,"Decay Electron Truth Eta");
0058 theTree->Branch("de_pt_truth",&de_pt_truth,"Decay Electron Truth Pt");
0059 theTree->Branch("de_p_truth",&de_p_truth,"Decay Electron Truth P");
0060 theTree->Branch("dp_phi_truth",&dp_phi_truth,"Decay Positron Truth Phi");
0061 theTree->Branch("dp_eta_truth",&dp_eta_truth,"Decay Positron Truth Eta");
0062 theTree->Branch("dp_pt_truth",&dp_pt_truth,"Decay Positron Truth Pt");
0063 theTree->Branch("dp_p_truth",&dp_p_truth,"Decay Positron Truth P");
0064 theTree->Branch("sp_phi_truth",&sp_phi_truth,"Scattered Proton Truth Phi");
0065 theTree->Branch("sp_eta_truth",&sp_eta_truth,"Scattered Proton Truth Eta");
0066 theTree->Branch("sp_pt_truth",&sp_pt_truth,"Scattered Proton Truth Pt");
0067 theTree->Branch("sp_p_truth",&sp_p_truth,"Scattered Proton Truth P");
0068 theTree->Branch("se_phi_truth",&se_phi_truth,"Scattered Electron Truth Phi");
0069 theTree->Branch("se_eta_truth",&se_eta_truth,"Scattered Electron Truth Eta");
0070 theTree->Branch("se_pt_truth",&se_pt_truth,"Scattered Electron Truth Pt");
0071 theTree->Branch("se_p_truth",&se_p_truth,"Scattered Electron Truth P");
0072 theTree->Branch("jpsi_phi_truth",&jpsi_phi_truth,"Jpsi Truth Phi");
0073 theTree->Branch("jpsi_eta_truth",&jpsi_eta_truth,"Jpsi Truth Eta");
0074 theTree->Branch("jpsi_pt_truth",&jpsi_pt_truth,"Jpsi Truth Pt");
0075 theTree->Branch("jpsi_p_truth",&jpsi_p_truth,"Jpsi Truth P");
0076 theTree->Branch("de_phi_reco",&de_phi_reco,"Decay Electron Reco Phi");
0077 theTree->Branch("de_eta_reco",&de_eta_reco,"Decay Electron Reco Eta");
0078 theTree->Branch("de_pt_reco",&de_pt_reco,"Decay Electron Reco Pt");
0079 theTree->Branch("de_p_reco",&de_p_reco,"Decay Electron Reco P");
0080 theTree->Branch("dp_phi_reco",&dp_phi_reco,"Decay Positron Reco Phi");
0081 theTree->Branch("dp_eta_reco",&dp_eta_reco,"Decay Positron Reco Eta");
0082 theTree->Branch("dp_pt_reco",&dp_pt_reco,"Decay Positron Reco Pt");
0083 theTree->Branch("dp_p_reco",&dp_p_reco,"Decay Positron Reco P");
0084 theTree->Branch("sp_phi_reco",&sp_phi_reco,"Scattered Proton Reco Phi");
0085 theTree->Branch("sp_eta_reco",&sp_eta_reco,"Scattered Proton Reco Eta");
0086 theTree->Branch("sp_pt_reco",&sp_pt_reco,"Scattered Proton Reco Pt");
0087 theTree->Branch("sp_p_reco",&sp_p_reco,"Scattered Proton Reco P");
0088 theTree->Branch("se_phi_reco",&se_phi_reco,"Scattered Electron Reco Phi");
0089 theTree->Branch("se_eta_reco",&se_eta_reco,"Scattered Electron Reco Eta");
0090 theTree->Branch("se_pt_reco",&se_pt_reco,"Scattered Electron Reco Pt");
0091 theTree->Branch("se_p_reco",&se_p_reco,"Scattered Electron Reco P");
0092 theTree->Branch("jpsi_phi_reco",&jpsi_phi_reco,"Jpsi Reco Phi");
0093 theTree->Branch("jpsi_eta_reco",&jpsi_eta_reco,"Jpsi Reco Eta");
0094 theTree->Branch("jpsi_pt_reco",&jpsi_pt_reco,"Jpsi Reco Pt");
0095 theTree->Branch("jpsi_p_reco",&jpsi_p_reco,"Jpsi Reco P");
0096 theTree->Branch("jpsi_px_truth",&jpsi_px_truth,"Jpsi Truth Px");
0097 theTree->Branch("jpsi_py_truth",&jpsi_py_truth,"Jpsi Truth Py");
0098 theTree->Branch("jpsi_pz_truth",&jpsi_pz_truth,"Jpsi Truth Pz");
0099 theTree->Branch("jpsi_px_reco",&jpsi_px_reco,"Jpsi Reco Px");
0100 theTree->Branch("jpsi_py_reco",&jpsi_py_reco,"Jpsi Reco Py");
0101 theTree->Branch("jpsi_pz_reco",&jpsi_pz_reco,"Jpsi Reco Pz");
0102 theTree->Branch("invariant_reco_decay",&invariant_reco_decay,"Invariant Mass Decay Particles");
0103 theTree->Branch("invariant_reco_scatter",&invariant_reco_scatter,"Invariant Mass Scattered Electron");
0104
0105
0106 TCut pcut("Smeared.particles.p<=50&&Smeared.particles.p>1");
0107
0108
0109 TCut ecut("-log(tan(Smeared.particles.theta/2))<=4&&-log(tan(Smeared.particles.theta/2))>=-4");
0110
0111
0112 TCut cut_scattered_lepton("particles.id==11&&particles.KS==2&&Smeared.particles.p<=50&&Smeared.particles.p>1");
0113
0114
0115 TCut cut_decay_electron("particles.id==11&&particles.KS==1&&Smeared.particles.p<=50&&Smeared.particles.p>1");
0116
0117
0118 TCut cut_decay_positron("particles.id==-11&&particles.KS==1&&Smeared.particles.p<=50&&Smeared.particles.p>1");
0119
0120
0121
0122
0123
0124 erhic::EventMilou *event = NULL;
0125 Smear::Event *eventS = NULL;
0126
0127 tree->SetBranchAddress("event",&event);
0128 tree_smeared->SetBranchAddress("eventS", &eventS);
0129
0130 unsigned max_event = tree_smeared->GetEntries();
0131 std::vector<float> particle_eta;
0132 particle_eta.push_back(0);
0133 particle_eta.push_back(0);
0134 particle_eta.push_back(0);
0135 particle_eta.push_back(0);
0136 particle_eta.push_back(0);
0137 std::vector<float> particle_theta;
0138 particle_theta.push_back(0);
0139 particle_theta.push_back(0);
0140 particle_theta.push_back(0);
0141 particle_theta.push_back(0);
0142 particle_theta.push_back(0);
0143 std::vector<float> particle_phi;
0144 particle_phi.push_back(0);
0145 particle_phi.push_back(0);
0146 particle_phi.push_back(0);
0147 particle_phi.push_back(0);
0148 particle_phi.push_back(0);
0149 std::vector<float> particle_pt;
0150 particle_pt.push_back(0);
0151 particle_pt.push_back(0);
0152 particle_pt.push_back(0);
0153 particle_pt.push_back(0);
0154 particle_pt.push_back(0);
0155 std::vector<float> particle_p;
0156 particle_p.push_back(0);
0157 particle_p.push_back(0);
0158 particle_p.push_back(0);
0159 particle_p.push_back(0);
0160 particle_p.push_back(0);
0161
0162 std::vector<float> true_particle_eta;
0163 true_particle_eta.push_back(0);
0164 true_particle_eta.push_back(0);
0165 true_particle_eta.push_back(0);
0166 true_particle_eta.push_back(0);
0167 true_particle_eta.push_back(0);
0168 std::vector<float> true_particle_theta;
0169 true_particle_theta.push_back(0);
0170 true_particle_theta.push_back(0);
0171 true_particle_theta.push_back(0);
0172 true_particle_theta.push_back(0);
0173 true_particle_theta.push_back(0);
0174 std::vector<float> true_particle_phi;
0175 true_particle_phi.push_back(0);
0176 true_particle_phi.push_back(0);
0177 true_particle_phi.push_back(0);
0178 true_particle_phi.push_back(0);
0179 true_particle_phi.push_back(0);
0180 std::vector<float> true_particle_p;
0181 true_particle_p.push_back(0);
0182 true_particle_p.push_back(0);
0183 true_particle_p.push_back(0);
0184 true_particle_p.push_back(0);
0185 true_particle_p.push_back(0);
0186 std::vector<float> true_particle_pt;
0187 true_particle_pt.push_back(0);
0188 true_particle_pt.push_back(0);
0189 true_particle_pt.push_back(0);
0190 true_particle_pt.push_back(0);
0191 true_particle_pt.push_back(0);
0192 for ( unsigned ievent = 0; ievent < max_event/1000; ievent++ )
0193 {
0194 if ( ievent%1000 == 0 )
0195 cout << "Processing event " << ievent << endl;
0196
0197
0198 tree->GetEntry(ievent);
0199 tree_smeared->GetEntry(ievent);
0200 cout << Q2fromEICTree << endl;
0201
0202 unsigned ntracks = eventS->GetNTracks();
0203
0204
0205
0206
0207 for ( unsigned itrack = 0; itrack < ntracks; itrack++ )
0208 {
0209
0210 Smear::ParticleMCS * sparticle = eventS->GetTrack( itrack );
0211 erhic::ParticleMC * tparticle = event->GetTrack(itrack);
0212 int type_of_particle = -1;
0213
0214
0215 if(tparticle->Id().Code()==11&&tparticle->GetStatus()==1)
0216 {
0217 type_of_particle = 0;
0218 }
0219
0220
0221 if(tparticle->Id().Code()==11&&tparticle->GetStatus()==2)
0222 {
0223 type_of_particle = 1;
0224 }
0225
0226
0227 if(tparticle->Id().Code()==-11&&tparticle->GetStatus()==1)
0228 {
0229 type_of_particle = 2;
0230 }
0231
0232
0233 if(tparticle->Id().Code()==2212&&tparticle->GetStatus()==1)
0234 {
0235 type_of_particle = 3;
0236 }
0237
0238 if(tparticle->Id().Code()==443)
0239 {
0240 type_of_particle = 4;
0241 }
0242
0243
0244 if(sparticle==NULL&&type_of_particle!=4)
0245 {
0246 particle_eta.at(type_of_particle)=NULL;
0247 particle_theta.at(type_of_particle)=NULL;
0248 particle_phi.at(type_of_particle)=NULL;
0249 particle_pt.at(type_of_particle)=NULL;
0250 particle_p.at(type_of_particle)=NULL;
0251
0252 true_particle_eta.at(type_of_particle)=tparticle->GetEta();
0253 true_particle_theta.at(type_of_particle)=tparticle->GetTheta();
0254 true_particle_phi.at(type_of_particle)=tparticle->GetPhi();
0255 true_particle_p.at(type_of_particle)=tparticle->GetP();
0256 true_particle_pt.at(type_of_particle)=tparticle->GetPt();
0257
0258
0259
0260
0261
0262 }
0263 else if(type_of_particle!=4)
0264 {
0265
0266 particle_eta.at(type_of_particle)=sparticle->GetEta();
0267 particle_theta.at(type_of_particle)=2*atan(exp(-sparticle->GetEta()));
0268 particle_phi.at(type_of_particle)=sparticle->GetPhi();
0269 particle_pt.at(type_of_particle)=sparticle->GetPt();
0270 particle_p.at(type_of_particle)=sparticle->GetP();
0271
0272 true_particle_eta.at(type_of_particle)=tparticle->GetEta();
0273 true_particle_theta.at(type_of_particle)=tparticle->GetTheta();
0274 true_particle_phi.at(type_of_particle)=tparticle->GetPhi();
0275 true_particle_p.at(type_of_particle)=tparticle->GetP();
0276 true_particle_pt.at(type_of_particle)=tparticle->GetPt();
0277 }
0278 else if(type_of_particle==4)
0279 {
0280 true_particle_eta.at(type_of_particle)=tparticle->GetEta();
0281 true_particle_theta.at(type_of_particle)=tparticle->GetTheta();
0282 true_particle_phi.at(type_of_particle)=tparticle->GetPhi();
0283 true_particle_p.at(type_of_particle)=tparticle->GetP();
0284 true_particle_pt.at(type_of_particle)=tparticle->GetPt();
0285 }
0286
0287 }
0288
0289
0290 bool can_we_reconstruct_3_final_state_leptons = true;
0291 if((particle_p.at(0)<1||particle_p.at(1)<1)||particle_p.at(2)<1)
0292 {
0293 can_we_reconstruct_3_final_state_leptons = false;
0294 }
0295
0296 bool can_we_reconstruct_scattered_lepton = true;
0297 if(particle_p.at(1)<1)
0298 {
0299 can_we_reconstruct_scattered_lepton = false;
0300 }
0301
0302 bool can_we_reconstruct_decay_particles = true;
0303 if(particle_p.at(0)<1||particle_p.at(2)<1)
0304 {
0305 can_we_reconstruct_decay_particles = false;
0306 }
0307
0308
0309
0310 if(can_we_reconstruct_3_final_state_leptons)
0311 {
0312 float m = sqrt(2*particle_pt.at(0)*particle_pt.at(2)*(cosh(particle_eta.at(0)-particle_eta.at(2))-cos(particle_phi.at(0)-particle_phi.at(2))));
0313 invariant_reco_decay=m;
0314 }
0315
0316
0317
0318 if(can_we_reconstruct_3_final_state_leptons)
0319 {
0320 float m = sqrt(2*particle_pt.at(1)*particle_pt.at(2)*(cosh(particle_eta.at(1)-particle_eta.at(2))-cos(particle_phi.at(1)-particle_phi.at(2))));
0321 invariant_reco_scatter=m;
0322 }
0323
0324
0325 if(can_we_reconstruct_3_final_state_leptons)
0326 {
0327 float sum_px = particle_pt.at(0)*cos(particle_phi.at(0))+particle_pt.at(2)*cos(particle_phi.at(2));
0328 float sum_py = particle_pt.at(0)*sin(particle_phi.at(0))+particle_pt.at(2)*sin(particle_phi.at(2));
0329 float sum_pz = particle_p.at(0)*cos(particle_theta.at(0))+particle_p.at(2)*cos(particle_theta.at(2));
0330 float sum_p = sqrt(sum_px*sum_px+sum_py*sum_py+sum_pz*sum_pz);
0331 jpsi_px_reco=sum_px;
0332 jpsi_py_reco=sum_py;
0333 jpsi_pz_reco=sum_pz;
0334
0335
0336 float sum_px_true = true_particle_pt.at(0)*cos(true_particle_phi.at(0))+true_particle_pt.at(2)*cos(true_particle_phi.at(2));
0337 float sum_py_true = true_particle_pt.at(0)*sin(true_particle_phi.at(0))+true_particle_pt.at(2)*sin(true_particle_phi.at(2));
0338 float sum_pz_true = true_particle_p.at(0)*cos(true_particle_theta.at(0))+true_particle_p.at(2)*cos(true_particle_theta.at(2));
0339 jpsi_px_truth=sum_px_true;
0340 jpsi_py_truth=sum_py_true;
0341 jpsi_pz_truth=sum_pz_true;
0342
0343
0344
0345
0346
0347
0348
0349
0350
0351 float jpsi_eta = -log(tan( (acos(jpsi_pz_reco/sum_p))/2 ) );
0352 float jpsi_phi = atan(jpsi_py_reco/jpsi_px_reco);
0353 if(jpsi_py_reco>0&&jpsi_px_reco<0)
0354 jpsi_phi=3.14159265+jpsi_phi;
0355 else if(jpsi_py_reco<0&&jpsi_px_reco<0)
0356 jpsi_phi=3.14159265+jpsi_phi;
0357 else if(jpsi_py_reco<0&&jpsi_px_reco>0)
0358 jpsi_phi=2*3.14159265+jpsi_phi;
0359
0360
0361
0362
0363 particle_p.at(4)=sum_p;
0364 particle_pt.at(4)=sum_p*sin(jpsi_phi);
0365 particle_eta.at(4)=jpsi_eta;
0366 particle_phi.at(4)=jpsi_phi;
0367
0368
0369
0370
0371
0372 t = 2*275*true_particle_p.at(3)*(1-cos(2*atan(exp(-true_particle_eta.at(3)))));
0373 }
0374
0375 de_phi_truth=true_particle_phi.at(0);
0376 de_eta_truth=true_particle_eta.at(0);
0377 de_pt_truth=true_particle_pt.at(0);
0378 de_p_truth=true_particle_p.at(0);
0379 dp_phi_truth=true_particle_phi.at(2);
0380 dp_eta_truth=true_particle_eta.at(2);
0381 dp_pt_truth=true_particle_pt.at(2);
0382 dp_p_truth=true_particle_p.at(2);
0383 sp_phi_truth=true_particle_phi.at(3);
0384 sp_eta_truth=true_particle_eta.at(3);
0385 sp_pt_truth=true_particle_pt.at(3);
0386 sp_p_truth=true_particle_p.at(3);
0387 se_phi_truth=true_particle_phi.at(1);
0388 se_eta_truth=true_particle_eta.at(1);
0389 se_pt_truth=true_particle_pt.at(1);
0390 se_p_truth=true_particle_p.at(1);
0391 jpsi_phi_truth=true_particle_phi.at(4);
0392 jpsi_eta_truth=true_particle_eta.at(4);
0393 jpsi_pt_truth=true_particle_pt.at(4);
0394 jpsi_p_truth=true_particle_p.at(4);
0395 de_phi_reco=particle_phi.at(0);
0396 de_eta_reco=particle_eta.at(0);
0397 de_pt_reco=particle_pt.at(0);
0398 de_p_reco=particle_p.at(0);
0399 dp_phi_reco=particle_phi.at(2);
0400 dp_eta_reco=particle_eta.at(2);
0401 dp_pt_reco=particle_pt.at(2);
0402 dp_p_reco=particle_p.at(2);
0403 sp_phi_reco=particle_phi.at(3);
0404 sp_eta_reco=particle_eta.at(3);
0405 sp_pt_reco=particle_pt.at(3);
0406 sp_p_reco=particle_p.at(3);
0407 se_phi_reco=particle_phi.at(1);
0408 se_eta_reco=particle_eta.at(1);
0409 se_pt_reco=particle_pt.at(1);
0410 se_p_reco=particle_p.at(1);
0411 jpsi_phi_reco=particle_phi.at(4);
0412 jpsi_eta_reco=particle_eta.at(4);
0413 jpsi_pt_reco=particle_pt.at(4);
0414 jpsi_p_reco=particle_p.at(4);
0415
0416
0417 theTree->Fill();
0418 }
0419
0420
0421
0422
0423
0424 myfile->cd();
0425 theTree->Write("Tree");
0426 myfile->Write();
0427 myfile->Close();
0428 return 0;
0429 }
0430
0431 float eta_to_theta(float eta)
0432 {
0433 return 2*atan(exp(-1*eta));
0434 }
0435
0436 float theta_to_eta(float theta)
0437 {
0438 return -log(tan(theta/2));
0439 }