Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MakeAmpToolsFlat.C
Go to the documentation of this file.
1 #define MakeAmpToolsFlat_cxx
2 #include "MakeAmpToolsFlat.h"
3 #include <TH2.h>
4 #include <TStyle.h>
5 #include <TCanvas.h>
6 
7 Int_t foption; // select which file to write
8 
9 
11 {
12  // This file was created in the following way
13  // root> TFile *f = TFile::Open("treeFlat_DSelector_Z2pi_trees.root")
14  // root> TTree *t = (TTree *) f->Get("pippimmisspb208_TreeFlat")
15  // root> t->MakeClass("MakeAmpToolsFlat_gen")
16  //
17 // In a ROOT session, you can do:
18 // root> .L MakeAmpToolsFlat.C
19 // root> MakeAmpToolsFlat t
20 // root> t.GetEntry(12); // Fill t data members with entry number 12
21 // root> t.Show(); // Show values of entry 12
22 // root> t.Show(16); // Read and show values of entry 16
23 // root> t.Loop(); // Loop on all entries
24 //
25 
26 // This is the loop skeleton where:
27 // jentry is the global entry number in the chain
28 // ientry is the entry number in the current Tree
29 // Note that the argument to GetEntry must be:
30 // jentry for TChain::GetEntry
31 // ientry for TTree::GetEntry and TBranch::GetEntry
32 //
33 // To read only selected branches, Insert statements like:
34 // METHOD1:
35 // fChain->SetBranchStatus("*",0); // disable all branches
36 // fChain->SetBranchStatus("branchname",1); // activate branchname
37 // METHOD2: replace line
38 // fChain->GetEntry(jentry); //read all branches
39 //by b_branchname->GetEntry(ientry); //read only this branch
40  if (fChain == 0) return;
41 
42 
43  // Define diagnostic histograms
44  TH1D *h1_twopimass_intime = new TH1D ("h1_twopimass_intime","Two pi mass intime (w=1)",200,0.2,0.8);
45  TH1D *h1_twopimass_outtime = new TH1D ("h1_twopimass_outtime","Two pi mass outtime (w=1)",200,0.2,0.8);
46  TH1D *h1_twopimass_total = new TH1D ("h1_twopimass_total","Two pi mass total (w=1)",200,0.2,0.8);
47  TH1D *h1_twopimass_signal = new TH1D ("h1_twopimass_signal","Two pi mass (weighted)",200,0.2,0.8);
48  TH1D *h1_twopimass_signalpos = new TH1D ("h1_twopimass_signalpos","Two pi mass (positive weighted)",200,0.2,0.8);
49 
50  // Three output Trees, depending on foption
51  // foption=1: First output has weights (positive and negative)
52  // foption=2: Second output has in-time only (weight=1)
53  // foption=3: Third output has out-time only (weight=negative 1/nbunches)
54  // foption=4: Third output has in-time only (weight=1 - out-time*(1/nbunches)/in-time)
55 
56  if (foption == 1){
57  outFile = new TFile("AmpToolsInputTree.root", "RECREATE");
58  }
59  else if (foption == 2) {
60  outFile = new TFile("AmpToolsInputTreeInTime.root", "RECREATE");
61  }
62  else if (foption == 3) {
63  outFile = new TFile("AmpToolsInputTreeOutTime.root", "RECREATE");
64  }
65  else if (foption == 4) {
66  outFile = new TFile("AmpToolsInputTreeInTimeW.root", "RECREATE");
67  }
68  else {
69  cout << "*** Loop, illegal foption=" << foption << endl;
70  }
71 
72  m_OutTree = new TTree("kin", "kin2");
73 
74  static size_t locNumFinalStateParticles = 3;
75 
76  m_OutTree->Branch("Weight", new float, "Weight/F");
77  m_OutTree->Branch("E_Beam", new float, "E_Beam/F");
78  m_OutTree->Branch("Px_Beam", new float, "Px_Beam/F");
79  m_OutTree->Branch("Py_Beam", new float, "Py_Beam/F");
80  m_OutTree->Branch("Pz_Beam", new float, "Pz_Beam/F");
81  m_OutTree->Branch("Target_Mass", new float, "Target_Mass/F");
82  m_OutTree->Branch("NumFinalState", new int, "NumFinalState/I");
83  m_OutTree->Branch("PID_FinalState", new int[locNumFinalStateParticles], "PID_FinalState[NumFinalState]/I");
84  m_OutTree->Branch("E_FinalState", new float[locNumFinalStateParticles], "E_FinalState[NumFinalState]/F");
85  m_OutTree->Branch("Px_FinalState", new float[locNumFinalStateParticles], "Px_FinalState[NumFinalState]/F");
86  m_OutTree->Branch("Py_FinalState", new float[locNumFinalStateParticles], "Py_FinalState[NumFinalState]/F");
87  m_OutTree->Branch("Pz_FinalState", new float[locNumFinalStateParticles], "Pz_FinalState[NumFinalState]/F");
88 
89  m_OutTree->SetBranchAddress("NumFinalState", &m_nPart);
90  m_nPart = 3;
91 
92  m_OutTree->SetBranchAddress("Target_Mass", &m_TargetMass);
93  m_TargetMass = 208*0.931494; // Pb mass in GeV.
94 
95  m_OutTree->SetBranchAddress("PID_FinalState", m_PID);
96  m_PID[0] = 8; m_PID[1] = 9; m_PID[2] = 111;
97 
98  m_OutTree->SetBranchAddress("E_FinalState", m_e);
99  m_OutTree->SetBranchAddress("Px_FinalState", m_px);
100  m_OutTree->SetBranchAddress("Py_FinalState", m_py);
101  m_OutTree->SetBranchAddress("Pz_FinalState", m_pz);
102  m_OutTree->SetBranchAddress("E_Beam", &m_eBeam);
103  m_OutTree->SetBranchAddress("Px_Beam", &m_pxBeam);
104  m_OutTree->SetBranchAddress("Py_Beam", &m_pyBeam);
105  m_OutTree->SetBranchAddress("Pz_Beam", &m_pzBeam);
106  m_OutTree->SetBranchAddress("Weight", &m_weight);
107 
108  // Process entries in Tree
109 
110  Long64_t nentries = fChain->GetEntriesFast();
111 
112  int NumIntime = 0;
113  int NumOuttime = 0;
114 
115  // loop over entries to count ratio of NumOuttime/NumIntime
116 
117  Long64_t nbytes = 0, nb = 0;
118  for (Long64_t jentry=0; jentry<nentries;jentry++) {
119  Long64_t ientry = LoadTree(jentry);
120  if (ientry < 0) break;
121  nb = fChain->GetEntry(jentry);
122  if (AccWeight > 0) NumIntime++;
123  if (AccWeight < 0) NumOuttime++;
124  }
125 
126  int Nbunches=4;
127  double signal_frac = (NumIntime - (1./Nbunches)*NumOuttime) / NumIntime;
128  cout << "Signal_frac=" << signal_frac << endl;
129 
130 
131  nbytes = 0, nb = 0;
132  for (Long64_t jentry=0; jentry<nentries;jentry++) {
133  Long64_t ientry = LoadTree(jentry);
134  if (ientry < 0) break;
135  nb = fChain->GetEntry(jentry);
136  /*nb = b_pip_p4_kin->GetEntry(jentry);
137  nbytes += nb;
138  nb = b_pim_p4_kin->GetEntry(jentry);
139  nbytes += nb;
140  nb = b_misspb_p4_kin->GetEntry(jentry);
141  nbytes += nb;
142  nb = b_beam_p4_kin->GetEntry(jentry);
143  nbytes += nb;
144  nb = b_AccWeight->GetEntry(jentry);*/
145  nbytes += nb;
146  // if (Cut(ientry) < 0) continue;
147 
148  /*cout<< "jentry=" << jentry << " px=" << beam_p4_kin->Px()<< " py=" << beam_p4_kin->Py()<< " pz=" << beam_p4_kin->Pz()<< " E=" << beam_p4_kin->E()<< endl;
149  cout<< "jentry=" << jentry << " weight=" << AccWeight << endl;
150  cout<< "jentry=" << jentry << " px=" << pip_p4_kin->Px()<< " py=" << pip_p4_kin->Py()<< " pz=" << pip_p4_kin->Pz()<< " E=" << pip_p4_kin->E()<< endl;
151  cout<< "jentry=" << jentry << " px=" << pim_p4_kin->Px()<< " py=" << pim_p4_kin->Py()<< " pz=" << pim_p4_kin->Pz()<< " E=" << pim_p4_kin->E()<< endl;
152  cout<< "jentry=" << jentry << " px=" << misspb_p4_kin->Px()<< " py=" << misspb_p4_kin->Py()<< " pz=" << misspb_p4_kin->Pz()<< " E=" << misspb_p4_kin->E()<< endl << endl;*/
153 
154  m_e[0] = pip_p4_kin->E();
155  m_px[0] = pip_p4_kin->Px();
156  m_py[0] = pip_p4_kin->Py();
157  m_pz[0] = pip_p4_kin->Pz();
158  m_e[1] = pim_p4_kin->E();
159  m_px[1] = pim_p4_kin->Px();
160  m_py[1] = pim_p4_kin->Py();
161  m_pz[1] = pim_p4_kin->Pz();
162  m_e[2] = misspb_p4_kin->E();
163  m_px[2] = misspb_p4_kin->Px();
164  m_py[2] = misspb_p4_kin->Py();
165  m_pz[2] = misspb_p4_kin->Pz();
166  m_eBeam = beam_p4_kin->E();
167  m_pxBeam = beam_p4_kin->Px();
168  m_pyBeam = beam_p4_kin->Py();
169  m_pzBeam = beam_p4_kin->Pz();
171  // m_weight = abs(AccWeight);
172 
173  TLorentzVector twopi;
174  twopi = *pip_p4_kin + *pim_p4_kin;
175  Double_t twopimass = twopi.M();
176 
177  if (foption == 1) {
178  if (AccWeight > 0) h1_twopimass_intime->Fill(twopimass);
179  if (AccWeight < 0) h1_twopimass_outtime->Fill(twopimass);
180  h1_twopimass_total->Fill(twopimass);
181  h1_twopimass_signal->Fill(twopimass,AccWeight);
182  h1_twopimass_signalpos->Fill(twopimass,abs(AccWeight));
183  m_OutTree->Fill();
184  }
185  else if (foption == 2) {
186  if (AccWeight > 0) m_OutTree->Fill(); // in time
187  }
188  else if (foption == 3) {
189  if (AccWeight < 0) {
190  m_weight = abs(AccWeight); // make all weights positve
191  m_OutTree->Fill(); // out of time
192  }
193  }
194  else if (foption == 4) {
195  if (AccWeight > 0) {
196  m_weight = signal_frac;
197  m_OutTree->Fill(); // in time, weight is fraction of signal
198  }
199  }
200  else {
201  cout << " *** Loop, illegal foption=" << foption << endl;
202  exit (1);
203  }
204  }
205 
206  // write out tree
207 
208  cout << "Completed loop: foption=" << foption << " nbytes =" << nbytes << " nentries=" << nentries << endl;
209  m_OutTree->Write();
210  outFile->Close();
211 
212  if (foption == 1) {
213  TFile *histFile = new TFile("AmpToolsDiagnosticHist.root", "RECREATE");
214  h1_twopimass_intime->Write();
215  h1_twopimass_outtime->Write();
216  h1_twopimass_total->Write();
217  h1_twopimass_signal->Write();
218  h1_twopimass_signalpos->Write();
219  histFile->Close();
220  }
221 
222 }
TLorentzVector * beam_p4_kin
TLorentzVector * pip_p4_kin
virtual Long64_t LoadTree(Long64_t entry)
virtual void Loop(Int_t foption)
Int_t foption
TLorentzVector * pim_p4_kin
TLorentzVector * misspb_p4_kin