18 dOutputFileName =
"DSelector_Z2pi_trees.root";
19 dOutputTreeFileName =
"tree_DSelector_Z2pi_trees.root";
20 dFlatTreeFileName =
"";
25 bool locInitializedPriorFlag = dInitializedFlag;
26 DSelector::Init(locTree);
28 if(locInitializedPriorFlag)
38 fMaxPion_dEdx =
new TF1(
"fMaxPion_dEdx",
"exp(-1.*[0]*x + [1]) + [2]", 0., 10.);
55 dAnalysisActions.push_back(
new DHistogramAction_ParticleID(dComboWrapper,
false));
56 dAnalysisActions.push_back(
new DHistogramAction_ParticleID(dComboWrapper,
true,
"KinFit"));
83 dHist_MissingMassSquared =
new TH1I(
"MissingMassSquared",
";Missing Mass Squared (GeV/c^{2})^{2}", 600, -0.24, 0.24);
84 dHist_BeamEnergy =
new TH1I(
"BeamEnergy",
";Beam Energy (GeV)", 600, 0.0, 12.0);
90 dHist_KinFitChiSq =
new TH1I(
"KinFitChiSq",
";Kinematic Fit #chi^{2}/NDF", 250, 0., 25.);
91 dHist_KinFitCL =
new TH1I(
"KinFitCL",
";Kinematic Fit Confidence Level", 100, 0., 1.);
93 dHist_M2pigen =
new TH1I(
"M2pigen",
";M_{#pi^{+}#pi^{-}} Gen (GeV/c^{2})", 400, 0.2, 0.6);
94 dHist_M2pikin =
new TH1I(
"M2pikin",
";M_{#pi^{+}#pi^{-}} Kin (GeV/c^{2})", 400, 0.2, 0.6);
95 dHist_M2pidiff =
new TH1I(
"M2pidiff",
";M_{#pi^{+}#pi^{-}} Kin - Gen (GeV/c^{2})", 400, -0.05, 0.05);
96 dHist_tgen =
new TH1I(
"tgen",
";|t| Gen (GeV/c)^{2}", 100, 0.0, 0.01);
97 dHist_tkin =
new TH1I(
"tkin",
";|t| Kin (GeV/c)^{2}", 100, 0.0, 0.01);
98 dHist_tdiff =
new TH1I(
"tdiff",
";|t| Kin - Gen (GeV/c)^{2}", 100, -0.01, 0.01);
99 dHist_tkin_tgen =
new TH2I(
"tkin_tgen",
"; |t| Gen ; |t| Kin (GeV/c)^{2}", 50, 0, 0.002, 50, 0, 0.002);
100 dHist_CosTheta_psi =
new TH2I(
"CosTheta_psi",
"; #psi; Cos#Theta", 90, -180., 180, 200, -1., 1.);
102 dHist_phigen_Phigen =
new TH2I(
"phigen_Phigen",
"; #Phi Gen; #phi Gen", 90, -180., 180, 90,-180,180);
103 dHist_phikin_Phikin =
new TH2I(
"phikin_Phikin",
"; #Phi Kin; #phi Kin", 90, -180., 180, 90,-180,180);
104 dHist_phimeas_phigen =
new TH2I(
"phimeas_phigen",
"; #phi Gen ; #phi Meas", 90, -180., 180, 90,-180,180);
105 dHist_phikin_phigen =
new TH2I(
"phikin_phigen",
"; #phi Gen ; #phi Kin", 90, -180., 180, 90,-180,180);
106 dHist_Phikin_Phigen =
new TH2I(
"Phikin_Phigen",
"; #Phi Gen ; #Phi Kin", 90, -180., 180, 90,-180,180);
107 dHist_Phimeas_Phigen =
new TH2I(
"Phimeas_Phigen",
"; #Phi Gen ; #Phi Meas", 90, -180., 180, 90,-180,180);
108 dHist_Delta_phi =
new TH2I(
"Delta_phi",
"; #phi ; #Delta #phi", 90, -180., 180, 90,-180,180);
109 dHist_Delta_Phi =
new TH2I(
"Delta_Phi",
"; #Phi ; #Delta #Phi", 90, -180., 180, 90,-180,180);
110 dHist_Delta_phimeas =
new TH2I(
"Delta_phimeas",
"; #phi ; #Delta #phi Meas", 90, -180., 180, 90,-180,180);
111 dHist_Delta_Phimeas =
new TH2I(
"Delta_Phimeas",
"; #Phi ; #Delta #Phi Meas", 90, -180., 180, 90,-180,180);
112 dHist_CosTheta =
new TH1I(
"CosTheta",
"; Cos#Theta", 100, -1., 1.);
113 dHist_CosThetadiff =
new TH1I(
"CosThetadiff",
"; Cos#Theta diff Kin-Gen", 50, -0.5, 0.5);
114 dHist_Phigen =
new TH1I(
"Phigen",
";Phigen (degrees)", 360,-180,180);
115 dHist_phigen =
new TH1I(
"phigen",
";phigen (degrees)", 360,-180,180);
116 dHist_Phikin =
new TH1I(
"Phikin",
";Phikin (degrees)", 360,-180,180);
117 dHist_phikin =
new TH1I(
"phikin",
";phikin (degrees)", 360,-180,180);
118 dHist_Phimeas =
new TH1I(
"Phimeas",
";Phimeas (degrees)", 360,-180,180);
119 dHist_phimeas =
new TH1I(
"phimeas",
";phimeas (degrees)", 360,-180,180);
120 dHist_psikin =
new TH1I(
"psikin",
";psi Kin (degrees)", 360,-180,180);
121 dHist_psigen =
new TH1I(
"psigen",
";psi Gen (degrees)", 360,-180,180);
122 dHist_Phidiff =
new TH1I(
"Phidiff",
";Phi Kin - Gen (degrees)", 100,-50,50);
123 dHist_phidiff =
new TH1I(
"phidiff",
";phi Kin - Gen (degrees)", 100,-50,50);
124 dHist_psidiff =
new TH1I(
"psidiff",
";psi Kin - Gen (degrees)", 100,-50,50);
126 dHist_pipDeltap =
new TH1I(
"pipDeltap",
"; #pi^{+}: Thrown p - KinFit p/Thrown p",100,-0.2,0.2);
127 dHist_pimDeltap =
new TH1I(
"pimDeltap",
"; #pi^{-}: Thrown p - KinFit p/ Thrown p",100,-0.2,0.2);
130 dHist_pimDeltap_Measured =
new TH1I(
"pimDeltap_Measured",
"; #pi^{-}: Thrown p - Measured p/ Thrown p",100,-0.2,0.2);
197 UInt_t locRunNumber = Get_RunNumber();
208 Reset_Actions_NewEvent();
218 set<Int_t> locUsedSoFar_BeamEnergy;
220 set<Int_t> locUsedSoFar_PiPlus, locUsedSoFar_PiMinus;
226 set<map<Particle_t, set<Int_t> > > locUsedSoFar_MissingMass, locUsedSoFar_2pi, locUsedSoFar_Angles;
248 double locEbeam_Thrown = 0;
249 if(dThrownBeam != NULL)
250 locEbeam_Thrown = dThrownBeam->Get_P4().E();
253 TLorentzVector locPb208P4_Thrown;
254 TLorentzVector locPiPlusP4_Thrown;
255 TLorentzVector locPiMinusP4_Thrown;
257 for(UInt_t loc_i = 0; loc_i < Get_NumThrown(); ++loc_i)
260 dThrownWrapper->Set_ArrayIndex(loc_i);
261 Particle_t thrown_pid = dThrownWrapper->Get_PID();
263 TLorentzVector locP4_Thrown = dThrownWrapper->Get_P4();
264 if (loc_i == 2) locPb208P4_Thrown = locP4_Thrown;
265 if (loc_i == 0) locPiPlusP4_Thrown = locP4_Thrown;
266 if (loc_i == 1) locPiMinusP4_Thrown = locP4_Thrown;
269 cout << endl <<
"Thrown" << endl;
270 cout <<
" dThrownBeam="; dThrownBeam->Get_P4().Print();
271 cout <<
" locPb208P4="; locPb208P4_Thrown.Print();
272 cout <<
" locPiPlusP4="; locPiPlusP4_Thrown.Print();
273 cout <<
" locPiMinusP4="; locPiMinusP4_Thrown.Print();
274 TLorentzVector loc2piP4_Thrown = locPiPlusP4_Thrown + locPiMinusP4_Thrown;
275 double tgen = (dThrownBeam->Get_P4() - locPiPlusP4_Thrown - locPiMinusP4_Thrown).M2();
283 for(UInt_t loc_i = 0; loc_i < Get_NumCombos(); ++loc_i)
286 dComboWrapper->Set_ComboIndex(loc_i);
289 if(dComboWrapper->Get_IsComboCut())
311 TLorentzVector locMissingP4 = locBeamP4;
312 locMissingP4 -= locPiPlusP4 + locPiMinusP4;
313 TLorentzVector loc2piP4 = locPiPlusP4 + locPiMinusP4;
315 cout <<
"Kin Fit" << endl;
316 cout <<
" locBeamP4="; locBeamP4.Print();
317 cout <<
" locMissingPb208P4="; locMissingPb208P4.Print();
318 cout <<
" locPiPlusP4="; locPiPlusP4.Print();
319 cout <<
" locPiMinusP4="; locPiMinusP4.Print();
320 cout <<
" locMissingP4="; locMissingP4.Print();
321 cout <<
" loc2piP4="; loc2piP4.Print();
324 cout <<
"KIN: MM2 =" << locMissingP4.M2() <<
" DeltaE=" << locBeamP4.E()-loc2piP4.E() <<
" GenDeltaE=" << locEbeam_Thrown-locBeamP4.E() << endl;
329 TLorentzVector locMissingPb208P4_Measured (0,0,0,193.750748);
331 TLorentzVector locPiPlusP4_Measured =
dPiPlusWrapper->Get_P4_Measured();
332 TLorentzVector locPiMinusP4_Measured =
dPiMinusWrapper->Get_P4_Measured();
334 TLorentzVector locMissingP4_Measured = locBeamP4_Measured;
335 locMissingP4_Measured -= locPiPlusP4_Measured + locPiMinusP4_Measured;
336 locMissingPb208P4_Measured += locMissingP4_Measured;
337 TLorentzVector loc2piP4_Measured = locPiPlusP4_Measured + locPiMinusP4_Measured;
339 cout <<
"Measured" << endl;
340 cout <<
" locBeamP4_Measured="; locBeamP4_Measured.Print();
341 cout <<
" locMissingPb208P4_Measured="; locMissingPb208P4_Measured.Print();
342 cout <<
" locPiPlusP4_Measured="; locPiPlusP4_Measured.Print();
343 cout <<
" locPiMinusP4_Measured="; locPiMinusP4_Measured.Print();
344 cout <<
" locMissingP4_Measured="; locMissingP4_Measured.Print();
345 cout <<
" loc2piP4_Measured="; loc2piP4_Measured.Print();
356 if(!Execute_Actions())
386 if(locPiPlus_dEdx_CDC >
fMaxPion_dEdx->Eval(locPiPlusP4.P()) || locPiMinus_dEdx_CDC >
fMaxPion_dEdx->Eval(locPiMinusP4.P())) {
387 dComboWrapper->Set_IsComboCut(
true);
391 cout <<
" Passed CDC dE/dX cut " << endl;
398 double locMissingMassSquared = locMissingP4_Measured.M2();
402 map<Particle_t, set<Int_t> > locUsedThisCombo_MissingMass;
403 locUsedThisCombo_MissingMass[
Unknown].insert(locBeamID);
404 locUsedThisCombo_MissingMass[
PiPlus].insert(locPiPlusTrackID);
405 locUsedThisCombo_MissingMass[
PiMinus].insert(locPiMinusTrackID);
409 dComboWrapper->Set_IsComboCut(
true);
413 cout <<
" Passed beam energy cut " << endl;
418 dComboWrapper->Set_IsComboCut(
true);
422 cout <<
" Passed Missing mass cut " << endl;
426 dHist_KinFitChiSq->Fill(dComboWrapper->Get_ChiSq_KinFit()/dComboWrapper->Get_NDF_KinFit());
427 dHist_KinFitCL->Fill(dComboWrapper->Get_ConfidenceLevel_KinFit());
428 if(dComboWrapper->Get_ConfidenceLevel_KinFit() <=
dMinKinFitCL) {
429 dComboWrapper->Set_IsComboCut(
true);
432 cout <<
" Passed kinematic Fit CL cut " << endl;
435 map<Particle_t, set<Int_t> > locUsedThisCombo_2piMass;
436 locUsedThisCombo_2piMass[
PiPlus].insert(locPiPlusTrackID);
437 locUsedThisCombo_2piMass[
PiMinus].insert(locPiMinusTrackID);
439 dComboWrapper->Set_IsComboCut(
true);
442 cout <<
" Passed 2pi mass cut " << endl;
445 if(locUsedSoFar_2pi.find(locUsedThisCombo_2piMass) == locUsedSoFar_2pi.end())
450 locUsedSoFar_2pi.insert(locUsedThisCombo_2piMass);
458 if(locUsedSoFar_BeamEnergy.find(locBeamID) == locUsedSoFar_BeamEnergy.end())
461 locUsedSoFar_BeamEnergy.insert(locBeamID);
465 if(locUsedSoFar_MissingMass.find(locUsedThisCombo_MissingMass) == locUsedSoFar_MissingMass.end())
469 locUsedSoFar_MissingMass.insert(locUsedThisCombo_MissingMass);
483 if(locUsedSoFar_PiPlus.find(locPiPlusTrackID) == locUsedSoFar_PiPlus.end())
486 dHist_pipDeltap->Fill(locPiPlusP4_Thrown.Vect().Mag() > 0? (locPiPlusP4_Thrown.Vect().Mag()-locPiPlusP4.Vect().Mag())/locPiPlusP4_Thrown.Vect().Mag() : 0);
487 dHist_pipDeltap_Measured->Fill(locPiPlusP4_Thrown.Vect().Mag() > 0? (locPiPlusP4_Thrown.Vect().Mag()-locPiPlusP4_Measured.Vect().Mag())/locPiPlusP4_Thrown.Vect().Mag() : 0);
488 locUsedSoFar_PiPlus.insert(locPiPlusTrackID);
491 if(locUsedSoFar_PiMinus.find(locPiMinusTrackID) == locUsedSoFar_PiMinus.end())
494 dHist_pimDeltap->Fill(locPiMinusP4_Thrown.Vect().Mag() > 0? (locPiMinusP4_Thrown.Vect().Mag()-locPiMinusP4.Vect().Mag())/locPiMinusP4_Thrown.Vect().Mag() : 0);
495 dHist_pimDeltap_Measured->Fill(locPiMinusP4_Thrown.Vect().Mag() > 0? (locPiMinusP4_Thrown.Vect().Mag()-locPiMinusP4_Measured.Vect().Mag())/locPiMinusP4_Thrown.Vect().Mag() : 0);
496 locUsedSoFar_PiMinus.insert(locPiMinusTrackID);
499 cout <<
"Ebeam Thrown P4.E=" << dThrownBeam->Get_P4().E() <<
" Kinfit P4.E=" << locBeamP4.E() <<
" P4 Measured =" << locBeamP4_Measured.E() << endl;
501 cout <<
"PiPlus Thrown P4.E=" << locPiPlusP4_Thrown.E() <<
" Kinfit P4.E=" << locPiPlusP4.E() <<
" P4 Measured =" << locPiPlusP4_Measured.E() << endl;
502 cout <<
"PiMinus Thrown P4.E=" << locPiMinusP4_Thrown.E() <<
" Kinfit P4.E=" << locPiMinusP4.E() <<
" P4 Measured =" << locPiMinusP4_Measured.E() << endl << endl;
505 if ( locBeamP4_Measured.Z() < 0 || locPiPlusP4_Measured.Z() < 0 || locPiMinusP4_Measured.Z() < 0 ||
506 locPiPlusP4.Z() < 0 || locPiMinusP4.Z() < 0 || (locPiPlusP4_Measured.Z() + locPiMinusP4_Measured.Z() + locMissingPb208P4_Measured_input.Z()) < 0 ||
507 (locPiPlusP4.Z() + locPiMinusP4.Z() + locMissingPb208P4.Z()) < 0 ) {
508 dComboWrapper->Set_IsComboCut(
true);
509 cout <<
"*** Failed Negative pz cut ***" << endl;
512 cout <<
" Passed Negative pz cut " << endl;
517 double tgen = (dThrownBeam->Get_P4() - loc2piP4_Thrown).M2();
518 TLorentzRotation resonanceBoost( -loc2piP4_Thrown.BoostVector() );
519 TLorentzVector beam_res = resonanceBoost * dThrownBeam->Get_P4();
520 TLorentzVector recoil_res = resonanceBoost * locPb208P4_Thrown;
521 TLorentzVector p1_res = resonanceBoost * locPiPlusP4_Thrown;
522 TLorentzVector p2_res = resonanceBoost * locPiMinusP4_Thrown;
525 TVector3 eps(cos(phipol),
sin(phipol), 0.0);
528 TVector3
y = (dThrownBeam->Get_P4().Vect().Unit().Cross(-locPb208P4_Thrown.Vect().Unit())).Unit();
531 TVector3 z = -1. * recoil_res.Vect().Unit();
532 TVector3
x = y.Cross(z).Unit();
533 TVector3 anglesgen( (p1_res.Vect()).Dot(x),
534 (p1_res.Vect()).Dot(y),
535 (p1_res.Vect()).Dot(z) );
537 double CosThetagen = anglesgen.CosTheta();
538 double phigen = anglesgen.Phi();
540 double Phigen = atan2(y.Dot(eps), dThrownBeam->Get_P4().Vect().Unit().Dot(eps.Cross(y)));
542 double psigen = Phigen - phigen;
543 if(psigen < -3.14159) psigen += 2*3.14159;
544 if(psigen > 3.14159) psigen -= 2*3.14159;
549 double tkin = (locBeamP4 - loc2piP4).M2();
550 TLorentzRotation resonanceBoost2( -loc2piP4.BoostVector() );
551 beam_res = resonanceBoost2 * locBeamP4;
552 recoil_res = resonanceBoost2 * locMissingPb208P4;
553 p1_res = resonanceBoost2 * locPiPlusP4;
554 p2_res = resonanceBoost2 * locPiMinusP4;
557 y = (locBeamP4.Vect().Unit().Cross(-locMissingPb208P4.Vect().Unit())).Unit();
560 z = -1. * recoil_res.Vect().Unit();
561 x = y.Cross(z).Unit();
562 TVector3 angleskin( (p1_res.Vect()).Dot(x),
563 (p1_res.Vect()).Dot(y),
564 (p1_res.Vect()).Dot(z) );
566 double CosThetakin = angleskin.CosTheta();
567 double phikin = angleskin.Phi();
569 double Phikin = atan2(y.Dot(eps), locBeamP4.Vect().Unit().Dot(eps.Cross(y)));
571 double psikin = Phikin - phikin;
572 if(psikin < -3.14159) psikin += 2*3.14159;
573 if(psikin > 3.14159) psikin -= 2*3.14159;
577 double tmeas = (locBeamP4_Measured - loc2piP4_Measured).M2();
578 TLorentzRotation resonanceBoost3( -loc2piP4_Measured.BoostVector() );
579 beam_res = resonanceBoost3 * locBeamP4_Measured;
580 recoil_res = resonanceBoost3 * locMissingPb208P4_Measured;
581 p1_res = resonanceBoost3 * locPiPlusP4_Measured;
582 p2_res = resonanceBoost3 * locPiMinusP4_Measured;
585 y = (locBeamP4_Measured.Vect().Unit().Cross(-locMissingPb208P4_Measured.Vect().Unit())).Unit();
588 z = -1. * recoil_res.Vect().Unit();
589 x = y.Cross(z).Unit();
590 TVector3 anglesmeas( (p1_res.Vect()).Dot(x),
591 (p1_res.Vect()).Dot(y),
592 (p1_res.Vect()).Dot(z) );
594 double CosThetameas = anglesmeas.CosTheta();
595 double phimeas = anglesmeas.Phi();
597 double Phimeas = atan2(y.Dot(eps), locBeamP4_Measured.Vect().Unit().Dot(eps.Cross(y)));
599 double psimeas = Phimeas - phimeas;
600 if(psimeas < -3.14159) psimeas += 2*3.14159;
601 if(psimeas > 3.14159) psimeas -= 2*3.14159;
603 double Delta_phi = phikin - phigen;
604 double Delta_Phi = Phikin - Phigen;
605 double Delta_phimeas = phimeas - phigen;
606 double Delta_Phimeas = Phimeas - Phigen;
608 map<Particle_t, set<Int_t> > locUsedThisCombo_Angles;
609 locUsedThisCombo_Angles[
Unknown].insert(locBeamID);
611 locUsedThisCombo_Angles[
PiPlus].insert(locPiPlusTrackID);
612 locUsedThisCombo_Angles[
PiMinus].insert(locPiMinusTrackID);
613 if(locUsedSoFar_Angles.find(locUsedThisCombo_Angles) == locUsedSoFar_Angles.end())
646 locUsedSoFar_Angles.insert(locUsedThisCombo_Angles);
673 Fill_NumCombosSurvivedHists();
722 Bool_t locIsEventCut =
true;
723 for(UInt_t loc_i = 0; loc_i < Get_NumCombos(); ++loc_i) {
725 dComboWrapper->Set_ComboIndex(loc_i);
727 if(dComboWrapper->Get_IsComboCut())
729 locIsEventCut =
false;
732 if(!locIsEventCut && dOutputTreeFileName !=
"")
751 DSelector::Finalize();
TH2I * dHist_phikin_phigen
TH1I * dHist_CosThetadiff
void Process(unsigned int &NEvents, unsigned int &NEvents_read)
TH1I * dHist_pipDeltap_Measured
DKinematicData * dMissingPb208Wrapper
TH2I * dHist_phikin_Phikin
TH2I * dHist_Phimeas_Phigen
TH1I * dHist_pimDeltap_Measured
DBeamParticle * dComboBeamWrapper
TH1I * dHist_piplusMomentumMeasured
TH1I * dHist_piminusMomentumMeasured
UInt_t dPreviousRunNumber
DChargedTrackHypothesis * dPiMinusWrapper
TH2I * dHist_phimeas_phigen
Double_t dMaxMissingMassSquared
TH2I * dHist_Delta_Phimeas
TH2I * dHist_CosThetakin_CosThetagen
TH2I * dHist_phigen_Phigen
Bool_t Process(Long64_t entry)
void Get_ComboWrappers(void)
DChargedTrackHypothesis * dPiPlusWrapper
TH2I * dHist_Phikin_Phigen
TH2I * dHist_Delta_phimeas
TH1I * dHist_MissingMassSquared
TH2I * dHist_CosTheta_psi
Double_t dMinMissingMassSquared