12 dOutputFileName =
"DSelector_Z2pi_trees.root";
13 dOutputTreeFileName =
"tree_DSelector_Z2pi_trees.root";
14 dFlatTreeFileName =
"treeFlat_DSelector_Z2pi_trees.root";
15 dFlatTreeName =
"pippimmisspb208_TreeFlat";
19 bool locInitializedPriorFlag = dInitializedFlag;
20 DSelector::Init(locTree);
22 if(locInitializedPriorFlag)
34 dAnalysisActions.push_back(
new DHistogramAction_ParticleID(dComboWrapper,
false));
35 dAnalysisActions.push_back(
new DHistogramAction_ParticleID(dComboWrapper,
true,
"KinFit"));
50 dAnalysisActions.push_back(
new DHistogramAction_BeamEnergy(dComboWrapper,
false));
63 dHist_MissingMassSquared =
new TH1I(
"MissingMassSquared",
";Missing Mass Squared (GeV/c^{2})^{2}", 600, -0.24, 0.24);
64 dHist_BeamEnergy =
new TH1I(
"BeamEnergy",
";Beam Energy (GeV)", 600, 0.0, 12.0);
70 dHist_KinFitChiSq =
new TH1I(
"KinFitChiSq",
";Kinematic Fit #chi^{2}/NDF", 250, 0., 25.);
71 dHist_KinFitCL =
new TH1I(
"KinFitCL",
";Kinematic Fit Confidence Level", 100, 0., 1.);
73 dHist_M2pigen =
new TH1I(
"M2pigen",
";M_{#pi^{+}#pi^{-}} Gen (GeV/c^{2})", 400, 0.2, 0.6);
74 dHist_M2pikin =
new TH1I(
"M2pikin",
";M_{#pi^{+}#pi^{-}} Kin (GeV/c^{2})", 400, 0.2, 0.6);
75 dHist_M2pidiff =
new TH1I(
"M2pidiff",
";M_{#pi^{+}#pi^{-}} Kin - Gen (GeV/c^{2})", 400, -0.05, 0.05);
76 dHist_tgen =
new TH1I(
"tgen",
";|t| Gen (GeV/c)^{2}", 100, 0.0, 0.01);
77 dHist_tkin =
new TH1I(
"tkin",
";|t| Kin (GeV/c)^{2}", 100, 0.0, 0.01);
78 dHist_tdiff =
new TH1I(
"tdiff",
";|t| Kin - Gen (GeV/c)^{2}", 100, -0.01, 0.01);
79 dHist_tkin_tgen =
new TH2I(
"tkin_tgen",
"; |t| Gen ; |t| Kin (GeV/c)^{2}", 50, 0, 0.002, 50, 0, 0.002);
80 dHist_CosTheta_psi =
new TH2I(
"CosTheta_psi",
"; #psi; Cos#Theta", 90, -180., 180, 200, -1., 1.);
82 dHist_phigen_Phigen =
new TH2I(
"phigen_Phigen",
"; #Phi Gen; #phi Gen", 90, -180., 180, 90,-180,180);
83 dHist_phikin_Phikin =
new TH2I(
"phikin_Phikin",
"; #Phi Kin; #phi Kin", 90, -180., 180, 90,-180,180);
84 dHist_phimeas_phigen =
new TH2I(
"phimeas_phigen",
"; #phi Gen ; #phi Meas", 90, -180., 180, 90,-180,180);
85 dHist_phikin_phigen =
new TH2I(
"phikin_phigen",
"; #phi Gen ; #phi Kin", 90, -180., 180, 90,-180,180);
86 dHist_Phikin_Phigen =
new TH2I(
"Phikin_Phigen",
"; #Phi Gen ; #Phi Kin", 90, -180., 180, 90,-180,180);
87 dHist_Phimeas_Phigen =
new TH2I(
"Phimeas_Phigen",
"; #Phi Gen ; #Phi Meas", 90, -180., 180, 90,-180,180);
88 dHist_Delta_phi =
new TH2I(
"Delta_phi",
"; #phi ; #Delta #phi", 90, -180., 180, 90,-180,180);
89 dHist_Delta_Phi =
new TH2I(
"Delta_Phi",
"; #Phi ; #Delta #Phi", 90, -180., 180, 90,-180,180);
90 dHist_Delta_phimeas =
new TH2I(
"Delta_phimeas",
"; #phi ; #Delta #phi Meas", 90, -180., 180, 90,-180,180);
91 dHist_Delta_Phimeas =
new TH2I(
"Delta_Phimeas",
"; #Phi ; #Delta #Phi Meas", 90, -180., 180, 90,-180,180);
92 dHist_CosTheta =
new TH1I(
"CosTheta",
"; Cos#Theta", 100, -1., 1.);
93 dHist_CosThetadiff =
new TH1I(
"CosThetadiff",
"; Cos#Theta diff Kin-Gen", 50, -0.5, 0.5);
94 dHist_Phigen =
new TH1I(
"Phigen",
";Phigen (degrees)", 360,-180,180);
95 dHist_phigen =
new TH1I(
"phigen",
";phigen (degrees)", 360,-180,180);
96 dHist_Phikin =
new TH1I(
"Phikin",
";Phikin (degrees)", 360,-180,180);
97 dHist_phikin =
new TH1I(
"phikin",
";phikin (degrees)", 360,-180,180);
98 dHist_Phimeas =
new TH1I(
"Phimeas",
";Phimeas (degrees)", 360,-180,180);
99 dHist_phimeas =
new TH1I(
"phimeas",
";phimeas (degrees)", 360,-180,180);
100 dHist_psikin =
new TH1I(
"psikin",
";psi Kin (degrees)", 360,-180,180);
101 dHist_psigen =
new TH1I(
"psigen",
";psi Gen (degrees)", 360,-180,180);
102 dHist_Phidiff =
new TH1I(
"Phidiff",
";Phi Kin - Gen (degrees)", 100,-50,50);
103 dHist_phidiff =
new TH1I(
"phidiff",
";phi Kin - Gen (degrees)", 100,-50,50);
104 dHist_psidiff =
new TH1I(
"psidiff",
";psi Kin - Gen (degrees)", 100,-50,50);
108 dHist_pipDeltap =
new TH1I(
"pipDeltap",
"; #pi^{+}: Thrown p - KinFit p/Thrown p",100,-0.2,0.2);
109 dHist_pimDeltap =
new TH1I(
"pimDeltap",
"; #pi^{-}: Thrown p - KinFit p/ Thrown p",100,-0.2,0.2);
112 dHist_pimDeltap_Measured =
new TH1I(
"pimDeltap_Measured",
"; #pi^{-}: Thrown p - Measured p/ Thrown p",100,-0.2,0.2);
119 fMaxPion_dEdx =
new TF1(
"fMaxPion_dEdx",
"exp(-1.*[0]*x + [1]) + [2]", 0., 10.);
156 dFlatTreeInterface->Create_Branch_Fundamental<Double_t>(
"AccWeight");
194 UInt_t locRunNumber = Get_RunNumber();
205 Reset_Actions_NewEvent();
215 set<Int_t> locUsedSoFar_BeamEnergy;
217 set<Int_t> locUsedSoFar_PiPlus, locUsedSoFar_PiMinus;
223 set<map<Particle_t, set<Int_t> > > locUsedSoFar_MissingMass, locUsedSoFar_2pi, locUsedSoFar_Angles;
244 double locEbeam_Thrown = 0;
245 if(dThrownBeam != NULL)
246 locEbeam_Thrown = dThrownBeam->Get_P4().E();
249 TLorentzVector locPb208P4_Thrown;
250 TLorentzVector locPiPlusP4_Thrown;
251 TLorentzVector locPiMinusP4_Thrown;
253 for(UInt_t loc_i = 0; loc_i < Get_NumThrown(); ++loc_i)
256 dThrownWrapper->Set_ArrayIndex(loc_i);
257 Particle_t thrown_pid = dThrownWrapper->Get_PID();
259 TLorentzVector locP4_Thrown = dThrownWrapper->Get_P4();
260 if (loc_i == 2) locPb208P4_Thrown = locP4_Thrown;
261 if (loc_i == 0) locPiPlusP4_Thrown = locP4_Thrown;
262 if (loc_i == 1) locPiMinusP4_Thrown = locP4_Thrown;
265 cout << endl <<
"Thrown" << endl;
266 cout <<
" dThrownBeam="; dThrownBeam->Get_P4().Print();
267 cout <<
" locPb208P4="; locPb208P4_Thrown.Print();
268 cout <<
" locPiPlusP4="; locPiPlusP4_Thrown.Print();
269 cout <<
" locPiMinusP4="; locPiMinusP4_Thrown.Print();
270 TLorentzVector loc2piP4_Thrown = locPiPlusP4_Thrown + locPiMinusP4_Thrown;
271 double tgen = (dThrownBeam->Get_P4() - locPiPlusP4_Thrown - locPiMinusP4_Thrown).M2();
277 cout <<
" Number of combos=" << Get_NumCombos() << endl;
278 for(UInt_t loc_i = 0; loc_i < Get_NumCombos(); ++loc_i)
281 dComboWrapper->Set_ComboIndex(loc_i);
284 if(dComboWrapper->Get_IsComboCut())
306 TLorentzVector locMissingP4 = locBeamP4;
307 locMissingP4 -= locPiPlusP4 + locPiMinusP4;
308 TLorentzVector loc2piP4 = locPiPlusP4 + locPiMinusP4;
311 cout << endl <<
"Kin Fit jcombo=" << loc_i << endl;
312 cout <<
" locBeamP4="; locBeamP4.Print();
313 cout <<
" locMissingPb208P4="; locMissingPb208P4.Print();
314 cout <<
" locPiPlusP4="; locPiPlusP4.Print();
315 cout <<
" locPiMinusP4="; locPiMinusP4.Print();
316 cout <<
" locMissingP4="; locMissingP4.Print();
317 cout <<
" loc2piP4="; loc2piP4.Print();
320 cout <<
"KIN: MM2 =" << locMissingP4.M2() <<
" DeltaE=" << locBeamP4.E()-loc2piP4.E() <<
" GenDeltaE=" << locEbeam_Thrown-locBeamP4.E() << endl;
325 TLorentzVector locMissingPb208P4_Measured (0,0,0,193.750748);
327 TLorentzVector locPiPlusP4_Measured =
dPiPlusWrapper->Get_P4_Measured();
328 TLorentzVector locPiMinusP4_Measured =
dPiMinusWrapper->Get_P4_Measured();
336 TLorentzVector locMissingP4_Measured = locBeamP4_Measured;
337 locMissingP4_Measured -= locPiPlusP4_Measured + locPiMinusP4_Measured;
338 locMissingPb208P4_Measured += locMissingP4_Measured;
339 TLorentzVector loc2piP4_Measured = locPiPlusP4_Measured + locPiMinusP4_Measured;
341 cout <<
"Measured" << endl;
342 cout <<
" locBeamP4_Measured="; locBeamP4_Measured.Print();
343 cout <<
" locMissingPb208P4_Measured="; locMissingPb208P4_Measured.Print();
344 cout <<
" locPiPlusP4_Measured="; locPiPlusP4_Measured.Print();
345 cout <<
" locPiMinusP4_Measured="; locPiMinusP4_Measured.Print();
346 cout <<
" locMissingP4_Measured="; locMissingP4_Measured.Print();
347 cout <<
" loc2piP4_Measured="; loc2piP4_Measured.Print();
353 if(!Execute_Actions())
387 if(locPiPlus_dEdx_CDC >
fMaxPion_dEdx->Eval(locPiPlusP4.P()) || locPiMinus_dEdx_CDC >
fMaxPion_dEdx->Eval(locPiMinusP4.P())) {
388 dComboWrapper->Set_IsComboCut(
true);
392 cout <<
" Passed CDC dE/dX cut " << endl;
397 double locMissingMassSquared = locMissingP4_Measured.M2();
401 map<Particle_t, set<Int_t> > locUsedThisCombo_MissingMass;
402 locUsedThisCombo_MissingMass[
Unknown].insert(locBeamID);
403 locUsedThisCombo_MissingMass[
PiPlus].insert(locPiPlusTrackID);
404 locUsedThisCombo_MissingMass[
PiMinus].insert(locPiMinusTrackID);
408 dComboWrapper->Set_IsComboCut(
true);
412 cout <<
" Passed beam energy cut " << endl;
417 dComboWrapper->Set_IsComboCut(
true);
421 cout <<
" Passed Missing mass cut " << endl;
424 dHist_KinFitChiSq->Fill(dComboWrapper->Get_ChiSq_KinFit(
"")/dComboWrapper->Get_NDF_KinFit(
""));
425 dHist_KinFitCL->Fill(dComboWrapper->Get_ConfidenceLevel_KinFit(
""));
426 if(dComboWrapper->Get_ConfidenceLevel_KinFit(
"") <=
dMinKinFitCL) {
427 dComboWrapper->Set_IsComboCut(
true);
430 cout <<
" Passed kinematic Fit CL cut " << endl;
433 map<Particle_t, set<Int_t> > locUsedThisCombo_2piMass;
434 locUsedThisCombo_2piMass[
PiPlus].insert(locPiPlusTrackID);
435 locUsedThisCombo_2piMass[
PiMinus].insert(locPiMinusTrackID);
437 dComboWrapper->Set_IsComboCut(
true);
440 cout <<
" Passed 2pi mass cut " << endl;
443 if(locUsedSoFar_2pi.find(locUsedThisCombo_2piMass) == locUsedSoFar_2pi.end())
448 locUsedSoFar_2pi.insert(locUsedThisCombo_2piMass);
456 if(locUsedSoFar_BeamEnergy.find(locBeamID) == locUsedSoFar_BeamEnergy.end())
459 locUsedSoFar_BeamEnergy.insert(locBeamID);
463 if(locUsedSoFar_MissingMass.find(locUsedThisCombo_MissingMass) == locUsedSoFar_MissingMass.end())
467 locUsedSoFar_MissingMass.insert(locUsedThisCombo_MissingMass);
481 if(locUsedSoFar_PiPlus.find(locPiPlusTrackID) == locUsedSoFar_PiPlus.end())
484 dHist_pipDeltap->Fill(locPiPlusP4_Thrown.Vect().Mag() > 0? (locPiPlusP4_Thrown.Vect().Mag()-locPiPlusP4.Vect().Mag())/locPiPlusP4_Thrown.Vect().Mag() : 0);
485 dHist_pipDeltap_Measured->Fill(locPiPlusP4_Thrown.Vect().Mag() > 0? (locPiPlusP4_Thrown.Vect().Mag()-locPiPlusP4_Measured.Vect().Mag())/locPiPlusP4_Thrown.Vect().Mag() : 0);
486 locUsedSoFar_PiPlus.insert(locPiPlusTrackID);
489 if(locUsedSoFar_PiMinus.find(locPiMinusTrackID) == locUsedSoFar_PiMinus.end())
492 dHist_pimDeltap->Fill(locPiMinusP4_Thrown.Vect().Mag() > 0? (locPiMinusP4_Thrown.Vect().Mag()-locPiMinusP4.Vect().Mag())/locPiMinusP4_Thrown.Vect().Mag() : 0);
493 dHist_pimDeltap_Measured->Fill(locPiMinusP4_Thrown.Vect().Mag() > 0? (locPiMinusP4_Thrown.Vect().Mag()-locPiMinusP4_Measured.Vect().Mag())/locPiMinusP4_Thrown.Vect().Mag() : 0);
494 locUsedSoFar_PiMinus.insert(locPiMinusTrackID);
497 cout <<
"Ebeam Thrown P4.E=" << dThrownBeam->Get_P4().E() <<
" Kinfit P4.E=" << locBeamP4.E() <<
" P4 Measured =" << locBeamP4_Measured.E() << endl;
499 cout <<
"PiPlus Thrown P4.E=" << locPiPlusP4_Thrown.E() <<
" Kinfit P4.E=" << locPiPlusP4.E() <<
" P4 Measured =" << locPiPlusP4_Measured.E() << endl;
500 cout <<
"PiMinus Thrown P4.E=" << locPiMinusP4_Thrown.E() <<
" Kinfit P4.E=" << locPiMinusP4.E() <<
" P4 Measured =" << locPiMinusP4_Measured.E() << endl;
503 if ( locBeamP4_Measured.Z() < 0 || locPiPlusP4_Measured.Z() < 0 || locPiMinusP4_Measured.Z() < 0 ||
504 locPiPlusP4.Z() < 0 || locPiMinusP4.Z() < 0 || locMissingPb208P4.Z() < 0 ) {
505 dComboWrapper->Set_IsComboCut(
true);
506 cout <<
"*** Failed Negative pz cut ***" << endl;
509 cout <<
" Passed Negative pz cut " << endl;
514 double tgen = (dThrownBeam->Get_P4() - loc2piP4_Thrown).M2();
515 TLorentzRotation resonanceBoost( -loc2piP4_Thrown.BoostVector() );
516 TLorentzVector beam_res = resonanceBoost * dThrownBeam->Get_P4();
517 TLorentzVector recoil_res = resonanceBoost * locPb208P4_Thrown;
518 TLorentzVector p1_res = resonanceBoost * locPiPlusP4_Thrown;
519 TLorentzVector p2_res = resonanceBoost * locPiMinusP4_Thrown;
522 TVector3 eps(cos(phipol),
sin(phipol), 0.0);
525 TVector3
y = (dThrownBeam->Get_P4().Vect().Unit().Cross(-locPb208P4_Thrown.Vect().Unit())).Unit();
528 TVector3 z = -1. * recoil_res.Vect().Unit();
529 TVector3
x = y.Cross(z).Unit();
530 TVector3 anglesgen( (p1_res.Vect()).Dot(x),
531 (p1_res.Vect()).Dot(y),
532 (p1_res.Vect()).Dot(z) );
534 double CosThetagen = anglesgen.CosTheta();
535 double phigen = anglesgen.Phi();
537 double Phigen = atan2(y.Dot(eps), dThrownBeam->Get_P4().Vect().Unit().Dot(eps.Cross(y)));
539 double psigen = Phigen - phigen;
540 if(psigen < -3.14159) psigen += 2*3.14159;
541 if(psigen > 3.14159) psigen -= 2*3.14159;
546 double tkin = (locBeamP4 - loc2piP4).M2();
547 TLorentzRotation resonanceBoost2( -loc2piP4.BoostVector() );
548 beam_res = resonanceBoost2 * locBeamP4;
549 recoil_res = resonanceBoost2 * locMissingPb208P4;
550 p1_res = resonanceBoost2 * locPiPlusP4;
551 p2_res = resonanceBoost2 * locPiMinusP4;
554 y = (locBeamP4.Vect().Unit().Cross(-locMissingPb208P4.Vect().Unit())).Unit();
557 z = -1. * recoil_res.Vect().Unit();
558 x = y.Cross(z).Unit();
559 TVector3 angleskin( (p1_res.Vect()).Dot(x),
560 (p1_res.Vect()).Dot(y),
561 (p1_res.Vect()).Dot(z) );
563 double CosThetakin = angleskin.CosTheta();
564 double phikin = angleskin.Phi();
566 double Phikin = atan2(y.Dot(eps), locBeamP4.Vect().Unit().Dot(eps.Cross(y)));
568 double psikin = Phikin - phikin;
569 if(psikin < -3.14159) psikin += 2*3.14159;
570 if(psikin > 3.14159) psikin -= 2*3.14159;
574 double tmeas = (locBeamP4_Measured - loc2piP4_Measured).M2();
575 TLorentzRotation resonanceBoost3( -loc2piP4_Measured.BoostVector() );
576 beam_res = resonanceBoost3 * locBeamP4_Measured;
577 recoil_res = resonanceBoost3 * locMissingPb208P4_Measured;
578 p1_res = resonanceBoost3 * locPiPlusP4_Measured;
579 p2_res = resonanceBoost3 * locPiMinusP4_Measured;
583 y = (locBeamP4_Measured.Vect().Unit().Cross(-locMissingPb208P4_Measured.Vect().Unit())).Unit();
586 z = -1. * recoil_res.Vect().Unit();
587 x = y.Cross(z).Unit();
588 TVector3 anglesmeas( (p1_res.Vect()).Dot(x),
589 (p1_res.Vect()).Dot(y),
590 (p1_res.Vect()).Dot(z) );
592 double CosThetameas = anglesmeas.CosTheta();
593 double phimeas = anglesmeas.Phi();
595 double Phimeas = atan2(y.Dot(eps), locBeamP4_Measured.Vect().Unit().Dot(eps.Cross(y)));
597 double psimeas = Phimeas - phimeas;
598 if(psimeas < -3.14159) psimeas += 2*3.14159;
599 if(psimeas > 3.14159) psimeas -= 2*3.14159;
601 double Delta_phi = phikin - phigen;
602 double Delta_Phi = Phikin - Phigen;
603 double Delta_phimeas = phimeas - phigen;
604 double Delta_Phimeas = Phimeas - Phigen;
607 map<Particle_t, set<Int_t> > locUsedThisCombo_Angles;
608 locUsedThisCombo_Angles[
Unknown].insert(locBeamID);
610 locUsedThisCombo_Angles[
PiPlus].insert(locPiPlusTrackID);
611 locUsedThisCombo_Angles[
PiMinus].insert(locPiMinusTrackID);
618 if(locUsedSoFar_Angles.find(locUsedThisCombo_Angles) == locUsedSoFar_Angles.end())
651 locUsedSoFar_Angles.insert(locUsedThisCombo_Angles);
661 double locRFTime = dComboWrapper->Get_RFTime_Measured();
665 double locBeamDeltaT = locBeam_X4_Measured.T() - (locRFTime + (locBeam_X4_Measured.Z() - dTargetCenter.Z())/29.9792458);
667 if(fabs(locBeamDeltaT) < 0.5*4.008) {
674 cout <<
" Tagger Accidentals: dTargetCenter=" << dTargetCenter.Z() <<
" locRFTime=" << locRFTime <<
" locBeamDeltaT=" << locBeamDeltaT <<
" AccWeight=" <<
AccWeight << endl;
675 cout <<
" locBeam_X4_Measured="; locBeam_X4_Measured.Print();
696 dFlatTreeInterface->Fill_Fundamental<Double_t>(
"AccWeight",
AccWeight);
699 dComboWrapper->Set_ComboIndex(loc_i);
704 Fill_NumCombosSurvivedHists();
753 Bool_t locIsEventCut =
true;
754 for(UInt_t loc_i = 0; loc_i < Get_NumCombos(); ++loc_i) {
756 dComboWrapper->Set_ComboIndex(loc_i);
758 if(dComboWrapper->Get_IsComboCut())
760 locIsEventCut =
false;
769 if (locBeamP4.Z() < 0 || locMissingPb208P4.Z() < 0 || locPiPlusP4.Z() < 0 || locPiMinusP4.Z() < 0 ) {
770 cout << endl <<
"Output: Negative Z momenta" << endl;
771 locBeamP4.Print(); locMissingPb208P4.Print(); locPiPlusP4.Print(); locPiMinusP4.Print();}
776 if(!locIsEventCut && dOutputTreeFileName !=
"")
797 DSelector::Finalize();
TH1I * dHist_piplusMomentumMeasured
TH2I * dHist_Delta_phimeas
TH2I * dHist_Phimeas_Phigen
void Process(unsigned int &NEvents, unsigned int &NEvents_read)
DChargedTrackHypothesis * dPiMinusWrapper
DBeamParticle * dComboBeamWrapper
TH2I * dHist_CosTheta_psi
Double_t dMaxMissingMassSquared
TH1I * dHist_TaggerAccidentals
Double_t dMinMissingMassSquared
DKinematicData * dMissingPb208Wrapper
void Get_ComboWrappers(void)
TH2I * dHist_CosThetakin_CosThetagen
TH1I * dHist_pipDeltap_Measured
TH2I * dHist_phimeas_phigen
TH1I * dHist_CosThetadiff
Bool_t Process(Long64_t entry)
TH2I * dHist_phikin_phigen
TH1I * dHist_pimDeltap_Measured
TH2I * dHist_phikin_Phikin
TH2I * dHist_Phikin_Phigen
TH2I * dHist_Delta_Phimeas
DChargedTrackHypothesis * dPiPlusWrapper
TH2I * dHist_phigen_Phigen
TH1I * dHist_MissingMassSquared
UInt_t dPreviousRunNumber
TH1I * dHist_piminusMomentumMeasured