]>
Commit | Line | Data |
---|---|---|
06f630bb | 1 | #define ana_sgl_cxx |
2 | #include "ana_sgl.h" | |
3 | #include <TH2.h> | |
4 | #include <TStyle.h> | |
5 | #include <TCanvas.h> | |
6 | #include <TMath.h> | |
7 | #include <vector> | |
8 | #include <deque> | |
9 | #include <cstdlib> | |
10 | #include <TRandom.h> | |
11 | ||
12 | using namespace std; | |
13 | ||
14 | using std::vector; | |
15 | ||
16 | class etrk; | |
17 | ||
18 | //Buffer for event mixing | |
19 | static const int NBUF=100; //depth of buffer | |
20 | static const int NMix=10; //# of events mixed (for +-) | |
21 | //static const int NMix=2; //# of events mixed (for +-) | |
22 | ||
23 | ||
24 | static const int NZBIN=10; | |
25 | static const int NCENT=10; | |
26 | int d_ibuf[NZBIN][NCENT]; | |
27 | vector<etrk> d_vep[NBUF][NZBIN][NCENT]; | |
28 | vector<etrk> d_vem[NBUF][NZBIN][NCENT]; | |
29 | ||
30 | static const unsigned int MAXPOOL=150; | |
31 | //static const unsigned int MAXPOOL=50; | |
32 | static const int MAX_TRY=3; | |
33 | deque<etrk> d_poolp[NZBIN][NCENT]; | |
34 | deque<etrk> d_poolm[NZBIN][NCENT]; | |
35 | ||
36 | void reshuffle_buffer(vector<etrk> &ve, | |
37 | deque<etrk> &pool){ | |
38 | //If there is not enough electron in the pool, give up | |
39 | unsigned int ne = ve.size(); | |
40 | unsigned int poolsize = pool.size(); | |
41 | if(poolsize < ne) { | |
42 | cout <<" pool size="<<poolsize<<" ne"<<ne<<endl; | |
43 | return; | |
44 | } | |
45 | for(unsigned int ie=0; ie < ne; ie++) { | |
46 | int j = rand()%poolsize; | |
47 | ve[ie] = pool[j]; | |
48 | } | |
49 | } | |
50 | ||
51 | void ana_sgl::Loop() | |
52 | { | |
53 | // In a ROOT session, you can do: | |
54 | // Root > .L ana_sgl.C | |
55 | // Root > ana_sgl t | |
56 | // Root > t.GetEntry(12); // Fill t data members with entry number 12 | |
57 | // Root > t.Show(); // Show values of entry 12 | |
58 | // Root > t.Show(16); // Read and show values of entry 16 | |
59 | // Root > t.Loop(); // Loop on all entries | |
60 | // | |
61 | ||
62 | // This is the loop skeleton where: | |
63 | // jentry is the global entry number in the chain | |
64 | // ientry is the entry number in the current Tree | |
65 | // Note that the argument to GetEntry must be: | |
66 | // jentry for TChain::GetEntry | |
67 | // ientry for TTree::GetEntry and TBranch::GetEntry | |
68 | // | |
69 | // To read only selected branches, Insert statements like: | |
70 | // METHOD1: | |
71 | // fChain->SetBranchStatus("*",0); // disable all branches | |
72 | // fChain->SetBranchStatus("branchname",1); // activate branchname | |
73 | // METHOD2: replace line | |
74 | // fChain->GetEntry(jentry); //read all branches | |
75 | //by b_branchname->GetEntry(ientry); //read only this branch | |
76 | if (fChain == 0) return; | |
77 | ||
78 | Long64_t nentries = fChain->GetEntriesFast(); | |
79 | ||
80 | Long64_t nbytes = 0, nb = 0; | |
81 | for (Long64_t jentry=0; jentry<nentries;jentry++) { | |
82 | Long64_t ientry = LoadTree(jentry); | |
83 | if (ientry < 0) break; | |
84 | nb = fChain->GetEntry(jentry); nbytes += nb; | |
85 | ana_event(ientry, jentry); | |
86 | // if (Cut(ientry) < 0) continue; | |
87 | } | |
88 | } | |
89 | ||
90 | //___________________________________________ | |
91 | void ana_sgl::ana_init(char *outname){ | |
92 | fout = new TFile(outname,"recreate"); | |
93 | char name[100]; | |
94 | ||
95 | hdedx_pt = new TH2D("hdedx_pt","dedx vs. pT",100,0, 10, 100, 0, 200); | |
96 | hdedx_tof_elec_pt = new TH2D("hdedx_tof_elec_pt","dedx vs. pT with tof veto",100,0, 10, 100, 0, 200); | |
97 | hdedx_tof_all_pt = new TH2D("hdedx_tof_all_pt","dedx vs. pT with tof veto",100,0, 10, 100, 0, 200); | |
98 | hdedx_tof_elec_emc_pt = new TH2D("hdedx_tof_elec_emc_pt","dedx vs. pT with tof veto and EMC",100,0, 10, 100, 0, 200); | |
99 | hdedx_tof_all_emc_pt = new TH2D("hdedx_tof_all_emc_pt","dedx vs. pT with tof veto and EMC",100,0, 10, 100, 0, 200); | |
100 | hdedx_emc_pt = new TH2D("hdedx_emc_pt","dedx vs. pT with EMC",100,0, 10, 100, 0, 200); | |
101 | ||
102 | ||
103 | hbetatof_pt = new TH2D("hbetatof_pt","beta vs. pT",100,0, 10, 100, 0, 1); | |
104 | hbetatof_tof_elec_pt = new TH2D("hbetatof_tof_elec_pt","beta vs. pT with tof cut",100,0, 10, 100, 0, 1); | |
105 | hbetatof_tof_all_pt = new TH2D("hbetatof_tof_all_pt","beta vs. pT with tof cut",100,0, 10, 100, 0, 1); | |
106 | hbetatof_tof_elec_emc_pt = new TH2D("hbetatof_tof_elec_emc_pt","beta vs. pT with tof cut",100,0, 10, 100, 0, 1); | |
107 | hbetatof_tof_all_emc_pt = new TH2D("hbetatof_tof_all_emc_pt","beta vs. pT with tof cut and emc E/p",100,0, 10, 100, 0, 1); | |
108 | hbetatof_emc_pt = new TH2D("hbetatof_emc_pt","dedx vs. pT with EMC",100,0, 10, 100, 0, 200); | |
109 | ||
110 | hCentrality = new TH1F("hCentrality","hCentrality",100,0,100); | |
111 | hV0AC = new TH2F("hV0AC","V0AC",5000,0,15000,5000,0,15000); | |
112 | hV0AC_Ntrk = new TH2F("hV0AC_Ntrk","V0AC vs. Ntrk",5000,0,30000,5000,0,15000); | |
113 | hV0AC_NaccTrcklts = new TH2F("hV0AC_NaccTrcklts","V0AC vs. Ntrklets",5000,0,30000,2500,0,5000); | |
114 | ||
115 | ||
116 | ||
117 | fBinWidth = 0.05; | |
118 | ||
119 | nHistos = (int)(10/fBinWidth); | |
120 | ||
121 | for(int i=0;i<nHistos;i++){ | |
122 | sprintf(name,"hdedx_p%d",i); | |
123 | hdedx[i] = new TH1D(name, name, 200,0,200); | |
124 | ||
125 | sprintf(name,"hdedx_tof_elec_p%d",i); | |
126 | hdedx_tof_elec[i] = new TH1D(name, name, 200,0,200); | |
127 | ||
128 | sprintf(name,"hdedx_tof_all_p%d",i); | |
129 | hdedx_tof_all[i] = new TH1D(name, name, 200,0,200); | |
130 | ||
131 | ||
132 | sprintf(name,"hdedx_tof_elec_emc_p%d",i); | |
133 | hdedx_tof_elec_emc[i] = new TH1D(name, name, 200,0,200); | |
134 | ||
135 | sprintf(name,"hdedx_tof_all_emc_p%d",i); | |
136 | hdedx_tof_all_emc[i] = new TH1D(name, name, 200,0,200); | |
137 | ||
138 | sprintf(name,"hdedx_emc_p%d",i); | |
139 | hdedx_emc[i] = new TH1D(name, name, 200,0,200); | |
140 | ||
141 | } | |
142 | ||
143 | ||
144 | int nbinx=400; | |
145 | float max_x=20; | |
146 | float min_x=0.2; | |
147 | float binw = (TMath::Log(max_x)-TMath::Log(min_x))/nbinx; | |
148 | double xbin[401]; | |
149 | for(int ii=0;ii<nbinx+1;ii++){ | |
150 | xbin[ii] = TMath::Exp(TMath::Log(min_x) + 0.5*binw+binw*ii); | |
151 | } | |
152 | ||
153 | ||
154 | //// my histo | |
155 | fEventStat = new TH1D("EventStat ","EventStat ",7,0,7); | |
156 | fEvent = new TH1D("Event","Number of Events", 60,0,60); | |
157 | fdEdXvsPt = new TH2D("dEdXvsPt","dE/dX vs. PT of TPC", nbinx, xbin, 2000,0,200); | |
158 | fdEdXnSigmaElecvsPt = new TH2D("fdEdXnSigmaElecvsPt","dE/dX normalized to electron vs. pT of TPC", | |
159 | nbinx, xbin, 2000, -10, 10); | |
160 | fTOFbetavsPt = new TH2D("fTOFbetavsPt","TOF beta vs. p", 400, 0, 20, 1200, 0, 1.2); | |
161 | fTOFnSigmaElecvsPt = new TH2D("fTOFnSigmaElecvsPt","TOF nsigma for electron", 400, 0, 20, 2000, -10, 10); | |
162 | ||
163 | fdEdXvsPtWithCut = new TH2D("dEdXvsPtWithCuts","dE/dX vs. PT of TPC", nbinx, xbin, 2000,0,200); | |
164 | ||
165 | ||
166 | fPtElec[0] = new TH1D("hPtElec0","elctron pt",250,0,5); | |
167 | fPtElec[1] = new TH1D("hPtElec1","elctron pt",250,0,5); | |
168 | ||
169 | ||
170 | fNelc_pos = new TH2D("hNelc_pos","# of electrons and positrons",40, -0.5, 39.5, 40, -0.5, 39.5); | |
171 | fNelc_all = new TH1D("hNelc_all","# of electrons and positrons",50, -0.5, 49.5); | |
172 | fNelc_all_pT = new TH2D("hNelc_all_pT","# of electrons and positrons in each pT bin",10, 0, 5, 50, -0.5, 49.5); | |
173 | ||
174 | for(int ic=0;ic<10;ic++){ | |
175 | sprintf(name,"hNelc_pos_cent%d",ic); | |
176 | fNelc_pos_cent[ic] = new TH2D(name, name, 40, -0.5, 39.5, 40, -0.5, 39.5); | |
177 | sprintf(name,"hNelc_all_cent%d",ic); | |
178 | fNelc_all_cent[ic] = new TH1D(name, name, 50, -0.5, 49.5); | |
179 | } | |
180 | ||
181 | ||
182 | //// tree branch | |
183 | d_tree = new TTree("d_tree","single tree"); | |
184 | d_tree->Branch("evt", &d_evt, "evt/I"); | |
185 | d_tree->Branch("cent", &d_cent, "cent/D"); | |
186 | d_tree->Branch("ntrk", &d_ntrk, "ntrk/D"); | |
187 | d_tree->Branch("xvprim", &d_xvprim, "xvprim/D"); | |
188 | d_tree->Branch("yvprim", &d_yvprim, "yvprim/D"); | |
189 | d_tree->Branch("zvprim", &d_zvprim, "zvprim/D"); | |
190 | d_tree->Branch("nacctrklets", &d_nacctrklets, "nacctrklets/D"); | |
191 | d_tree->Branch("xres", &d_xres, "xres/D"); | |
192 | d_tree->Branch("yres", &d_yres, "yres/D"); | |
193 | d_tree->Branch("zres", &d_zres, "zres/D"); | |
194 | d_tree->Branch("nelc", &d_nelc, "nelc/I"); | |
195 | d_tree->Branch("px", d_px, "px[nelc]/D"); | |
196 | d_tree->Branch("py", d_py, "py[nelc]/D"); | |
197 | d_tree->Branch("pz", d_pz, "pz[nelc]/D"); | |
198 | d_tree->Branch("p", d_p, "p[nelc]/D"); | |
199 | d_tree->Branch("pt", d_pt, "pt[nelc]/D"); | |
200 | d_tree->Branch("xv", d_xv, "xv[nelc]/D"); | |
201 | d_tree->Branch("yv", d_yv, "yv[nelc]/D"); | |
202 | d_tree->Branch("zv", d_zv, "zv[nelc]/D"); | |
203 | d_tree->Branch("phi", d_phi, "phi[nelc]/D"); | |
204 | d_tree->Branch("theta", d_theta, "theta[nelc]/D"); | |
205 | d_tree->Branch("eta", d_eta, "eta[nelc]/D"); | |
206 | d_tree->Branch("c", d_c, "c[nelc]/D"); | |
207 | d_tree->Branch("nclusITS", d_nclusITS, "nclusITS[nelc]/D"); | |
208 | d_tree->Branch("nclusTPC", d_nclusTPC, "nclusTPC[nelc]/D"); | |
209 | d_tree->Branch("nclusTPCiter", d_nclusTPCiter, "nclusTPCiter[nelc]/D"); | |
210 | d_tree->Branch("nfclusTPC", d_nfclusTPC, "nfclusTPC[nelc]/D"); | |
211 | d_tree->Branch("nfclusTPCr", d_nfclusTPCr, "nfclusTPCr[nelc]/D"); | |
212 | d_tree->Branch("nfclusTPCrFrac", d_nfclusTPCrFrac, "nfclusTPCrFrac[nelc]/D"); | |
213 | d_tree->Branch("TPCsignalN", d_TPCsignalN, "TPCsignalN[nelc]/D"); | |
214 | d_tree->Branch("TPCsignalNfrac", d_TPCsignalNfrac, "TPCsignalNfrac[nelc]/D"); | |
215 | d_tree->Branch("TPCchi2cl", d_TPCchi2cl, "TPCchi2cl[nelc]/D"); | |
216 | d_tree->Branch("trkstat", d_trkstat, "trkstat[nelc]/D"); | |
217 | d_tree->Branch("nclsTRD", d_nclsTRD, "nclsTRD[nelc]/D"); | |
218 | d_tree->Branch("TRDntracklets", d_TRDntracklets, "TRDntracklets[nelc]/D"); | |
219 | d_tree->Branch("TRDpidquality", d_TRDpidquality, "TRDpidquality[nelc]/D"); | |
220 | d_tree->Branch("TRDprobEle", d_TRDprobEle, "TRDprobEle[nelc]/D"); | |
221 | d_tree->Branch("TRDprobPio", d_TRDprobPio, "TRDprobPio[nelc]/D"); | |
222 | d_tree->Branch("impactXY", d_impactXY, "impactXY[nelc]/D"); | |
223 | d_tree->Branch("impactZ", d_impactZ, "impactZ[nelc]/D"); | |
224 | d_tree->Branch("tracklength", d_tracklength, "tracklength[nelc]/D"); | |
225 | d_tree->Branch("ITSsignal", d_ITSsignal, "ITSsignal[nelc]/D"); | |
226 | d_tree->Branch("ITSnsigmaEle", d_ITSnsigmaEle, "ITSnsigmaEle[nelc]/D"); | |
227 | d_tree->Branch("ITSnsigmaPio", d_ITSnsigmaPio, "ITSnsigmaPio[nelc]/D"); | |
228 | d_tree->Branch("ITSnsigmaMuo", d_ITSnsigmaMuo, "ITSnsigmaMuo[nelc]/D"); | |
229 | d_tree->Branch("ITSnsigmaKao", d_ITSnsigmaKao, "ITSnsigmaKao[nelc]/D"); | |
230 | d_tree->Branch("ITSnsigmaPro", d_ITSnsigmaPro, "ITSnsigmaPro[nelc]/D"); | |
231 | d_tree->Branch("PIn", d_PIn, "PIn[nelc]/D"); | |
232 | d_tree->Branch("TPCsignal", d_TPCsignal, "TPCsignal[nelc]/D"); | |
233 | d_tree->Branch("TOFsignal", d_TOFsignal, "TOFsignal[nelc]/D"); | |
234 | d_tree->Branch("TOFbeta", d_TOFbeta, "TOFbeta[nelc]/D"); | |
235 | d_tree->Branch("TPCnSigmaEle", d_TPCnSigmaEle, "TPCnSigmaEle[nelc]/D"); | |
236 | d_tree->Branch("TPCnSigmaPio", d_TPCnSigmaPio, "TPCnSigmaPio[nelc]/D"); | |
237 | d_tree->Branch("TPCnSigmaMuo", d_TPCnSigmaMuo, "TPCnSigmaMuo[nelc]/D"); | |
238 | d_tree->Branch("TPCnSigmaKao", d_TPCnSigmaKao, "TPCnSigmaKao[nelc]/D"); | |
239 | d_tree->Branch("TPCnSigmaPro", d_TPCnSigmaPro, "TPCnSigmaPro[nelc]/D"); | |
240 | d_tree->Branch("TOFnSigmaEle", d_TOFnSigmaEle, "TOFnSigmaEle[nelc]/D"); | |
241 | d_tree->Branch("TOFnSigmaPio", d_TOFnSigmaPio, "TOFnSigmaPio[nelc]/D"); | |
242 | d_tree->Branch("TOFnSigmaMuo", d_TOFnSigmaMuo, "TOFnSigmaMuo[nelc]/D"); | |
243 | d_tree->Branch("TOFnSigmaKao", d_TOFnSigmaKao, "TOFnSigmaKao[nelc]/D"); | |
244 | d_tree->Branch("TOFnSigmaPro", d_TOFnSigmaPro, "TOFnSigmaPro[nelc]/D"); | |
245 | d_tree->Branch("E", d_E, "E[nelc]/D"); | |
246 | d_tree->Branch("dphi", d_phi, "dphi[nelc]/D"); | |
247 | d_tree->Branch("deta", d_eta, "deta[nelc]/D"); | |
248 | ||
249 | d_tree->Branch("chi2ndf", d_chi2ndf, "chi2ndf[nelc]/D"); | |
250 | ||
251 | cout<<"aaaaaaaaaaa"<<endl; | |
252 | d_ntpair = new TTree("ntpair","pair"); | |
253 | d_ntpair->Branch("run", &d_run, "run/D"); | |
254 | d_ntpair->Branch("event", &d_event, "event/I"); | |
255 | d_ntpair->Branch("centrality",&d_centrality,"centrality/D"); | |
256 | d_ntpair->Branch("prim_xv",&d_prim_xv,"prim_xv/D"); | |
257 | d_ntpair->Branch("prim_yv",&d_prim_yv,"prim_yv/D"); | |
258 | d_ntpair->Branch("prim_zv",&d_prim_zv,"prim_zv/D"); | |
259 | d_ntpair->Branch("mass",&d_mass,"mass/D"); | |
260 | d_ntpair->Branch("pxpair",&d_pxpair,"pxpair/D"); | |
261 | d_ntpair->Branch("pypair",&d_pypair,"pypair/D"); | |
262 | d_ntpair->Branch("pzpair",&d_pzpair,"pzpair/D"); | |
263 | d_ntpair->Branch("ptpair",&d_ptpair,"ptpair/D"); | |
264 | d_ntpair->Branch("epair",&d_epair,"epair/D"); | |
265 | d_ntpair->Branch("eta",&d_etapair,"eta/D"); | |
266 | d_ntpair->Branch("phi",&d_phipair,"phi/D"); | |
267 | d_ntpair->Branch("cos",&d_cos,"cos/D"); | |
268 | d_ntpair->Branch("phiv",&d_phiv,"phiv/D"); | |
269 | d_ntpair->Branch("psi",&d_psi,"psi/D"); | |
270 | d_ntpair->Branch("pairtype",&d_pairtype,"pairtype/I"); | |
271 | d_ntpair->Branch("cent1",&d_cent1,"cent1/D"); | |
272 | d_ntpair->Branch("xv1",&d_xv1,"xv1/D"); | |
273 | d_ntpair->Branch("yv1",&d_yv1,"yv1/D"); | |
274 | d_ntpair->Branch("zv1",&d_zv1,"zv1/D"); | |
275 | d_ntpair->Branch("px1",&d_px1,"px1/D"); | |
276 | d_ntpair->Branch("py1",&d_py1,"py1/D"); | |
277 | d_ntpair->Branch("pz1",&d_pz1,"pz1/D"); | |
278 | d_ntpair->Branch("pt1",&d_pt1,"pt1/D"); | |
279 | d_ntpair->Branch("eta1",&d_eta1,"eta1/D"); | |
280 | d_ntpair->Branch("phi1",&d_phi1,"phi1/D"); | |
281 | d_ntpair->Branch("theta1",&d_theta1,"theta1/D"); | |
282 | d_ntpair->Branch("tpc1",&d_tpc1,"tpc1/D"); | |
283 | d_ntpair->Branch("ntpc_ele1",&d_ntpc_ele1,"ntpc_ele1/D"); | |
284 | d_ntpair->Branch("ntpc_pio1",&d_ntpc_pio1,"ntpc_pio1/D"); | |
285 | d_ntpair->Branch("ntpc_kao1",&d_ntpc_kao1,"ntpc_kao1/D"); | |
286 | d_ntpair->Branch("ntpc_pro1",&d_ntpc_pro1,"ntpc_pro1/D"); | |
287 | d_ntpair->Branch("beta1",&d_beta1,"beta1/D"); | |
288 | d_ntpair->Branch("ntof_ele1",&d_ntof_ele1,"ntof_ele1/D"); | |
289 | d_ntpair->Branch("ntof_pio1",&d_ntof_pio1,"ntof_pio1/D"); | |
290 | d_ntpair->Branch("ntof_kao1",&d_ntof_kao1,"ntof_kao1/D"); | |
291 | d_ntpair->Branch("ntof_pro1",&d_ntof_pro1,"ntof_pro1/D"); | |
292 | d_ntpair->Branch("its1",&d_its1,"its1/D"); | |
293 | d_ntpair->Branch("nits1",&d_nits1,"nits1/D"); | |
294 | d_ntpair->Branch("ntpc1",&d_ntpc1,"ntpc1/D"); | |
295 | d_ntpair->Branch("e1",&d_e1,"e1/D"); | |
296 | d_ntpair->Branch("dphi1",&d_dphi1,"dphi1/D"); | |
297 | d_ntpair->Branch("deta1",&d_deta1,"deta1/D"); | |
298 | d_ntpair->Branch("dcaxy1",&d_dcaxy1,"dcaxy1/D"); | |
299 | d_ntpair->Branch("dcaz1",&d_dcaz1,"dcaz1/D"); | |
300 | d_ntpair->Branch("conv1",&d_conv1,"conv1/I"); | |
301 | ||
302 | ||
303 | ||
304 | ||
305 | d_ntpair->Branch("cent2",&d_cent2,"cent2/D"); | |
306 | d_ntpair->Branch("xv2",&d_xv2,"xv2/D"); | |
307 | d_ntpair->Branch("yv2",&d_yv2,"yv2/D"); | |
308 | d_ntpair->Branch("zv2",&d_zv2,"zv2/D"); | |
309 | d_ntpair->Branch("px2",&d_px2,"px2/D"); | |
310 | d_ntpair->Branch("py2",&d_py2,"py2/D"); | |
311 | d_ntpair->Branch("pz2",&d_pz2,"pz2/D"); | |
312 | d_ntpair->Branch("pt2",&d_pt2,"pt2/D"); | |
313 | d_ntpair->Branch("eta2",&d_eta2,"eta2/D"); | |
314 | d_ntpair->Branch("phi2",&d_phi2,"phi2/D"); | |
315 | d_ntpair->Branch("theta2",&d_theta2,"theta2/D"); | |
316 | d_ntpair->Branch("tpc2",&d_tpc2,"tpc2/D"); | |
317 | d_ntpair->Branch("ntpc_ele2",&d_ntpc_ele2,"ntpc_ele2/D"); | |
318 | d_ntpair->Branch("ntpc_pio2",&d_ntpc_pio2,"ntpc_pio2/D"); | |
319 | d_ntpair->Branch("ntpc_kao2",&d_ntpc_kao2,"ntpc_kao2/D"); | |
320 | d_ntpair->Branch("ntpc_pro2",&d_ntpc_pro2,"ntpc_pro2/D"); | |
321 | d_ntpair->Branch("beta2",&d_beta2,"beta2/D"); | |
322 | d_ntpair->Branch("ntof_ele2",&d_ntof_ele2,"ntof_ele2/D"); | |
323 | d_ntpair->Branch("ntof_pio2",&d_ntof_pio2,"ntof_pio2/D"); | |
324 | d_ntpair->Branch("ntof_kao2",&d_ntof_kao2,"ntof_kao2/D"); | |
325 | d_ntpair->Branch("ntof_pro2",&d_ntof_pro2,"ntof_pro2/D"); | |
326 | d_ntpair->Branch("its2",&d_its2,"its2/D"); | |
327 | d_ntpair->Branch("nits2",&d_nits2,"nits2/D"); | |
328 | d_ntpair->Branch("ntpc2",&d_ntpc2,"ntpc2/D"); | |
329 | d_ntpair->Branch("e2",&d_e2,"e2/D"); | |
330 | d_ntpair->Branch("dphi2",&d_dphi2,"dphi2/D"); | |
331 | d_ntpair->Branch("deta2",&d_deta2,"deta2/D"); | |
332 | d_ntpair->Branch("dcaxy2",&d_dcaxy2,"dcaxy2/D"); | |
333 | d_ntpair->Branch("dcaz2",&d_dcaz2,"dcaz2/D"); | |
334 | d_ntpair->Branch("conv2",&d_conv2,"conv2/I"); | |
335 | ||
336 | cout<<"aaaaaaaaaaa"<<endl; | |
337 | for(int i=0;i<7;i++){ | |
338 | for(int j=0;j<11;j++){ | |
339 | sprintf(name,"hmasspt_cent%d_pair%d",j,i); | |
340 | hmasspt[i][j] = new TH2D(name, name, 500, 0, 5, 500, 0, 5); | |
341 | sprintf(name,"hmasspt_weight_cent%d_pair%d",j,i); | |
342 | hmasspt_weight[i][j] = new TH2D(name, name, 500, 0, 5, 500, 0, 5); | |
343 | } | |
344 | } | |
345 | ||
346 | vep.clear(); | |
347 | vem.clear(); | |
348 | vep_tmp.clear(); | |
349 | vem_tmp.clear(); | |
350 | ||
351 | d_evt = 0; | |
352 | d_event = 0; | |
353 | ||
354 | simflag = false; | |
355 | d_conv_flag = false; | |
356 | ||
357 | for(int i=0;i<10;i++){ | |
358 | nelec_pos[i] = 0; | |
359 | } | |
360 | ||
361 | magnetic_field_mm = true; | |
362 | ||
363 | d_flag_tof_cut = false; | |
364 | d_flag_emc_cut =false; | |
365 | d_flag_phiv = false; | |
366 | d_tpc_dedx_low = 75; | |
367 | d_tpc_dedx_high = 90; | |
368 | d_tof_low = -3; | |
369 | d_tof_high = 3; | |
370 | d_emc_low = 0.7; | |
371 | d_emc_high = 1.3; | |
372 | d_phiv_cut = 0.6; | |
373 | ||
374 | d_flag_kaon_veto = false; | |
375 | d_flag_proton_veto = false; | |
376 | d_dedx_kaon_veto_low = -2; | |
377 | d_dedx_kaon_veto_high = 2; | |
378 | d_dedx_proton_veto_low = -2; | |
379 | d_dedx_proton_veto_high = 2; | |
380 | ||
381 | ||
382 | cout<<"ana_init end:"<<endl; | |
383 | } | |
384 | ||
385 | //___________________________________________ | |
386 | void ana_sgl::ana_end(void){ | |
387 | ||
388 | fout->cd(); | |
389 | /* | |
390 | hdedx_pt->Write(); | |
391 | hdedx_tof_elec_pt->Write(); | |
392 | hdedx_tof_all_pt->Write(); | |
393 | hdedx_tof_elec_emc_pt->Write(); | |
394 | hdedx_tof_all_emc_pt->Write(); | |
395 | hdedx_emc_pt->Write(); | |
396 | ||
397 | hbetatof_pt->Write(); | |
398 | hbetatof_tof_elec_pt->Write(); | |
399 | hbetatof_tof_all_pt->Write(); | |
400 | hbetatof_tof_elec_emc_pt->Write(); | |
401 | hbetatof_tof_all_emc_pt->Write(); | |
402 | hbetatof_emc_pt->Write(); | |
403 | ||
404 | for(int i=0;i<nHistos;i++){ | |
405 | hdedx[i]->Write(); | |
406 | hdedx_tof_elec[i]->Write(); | |
407 | hdedx_tof_all[i]->Write(); | |
408 | hdedx_tof_elec_emc[i]->Write(); | |
409 | hdedx_tof_all_emc[i]->Write(); | |
410 | hdedx_emc[i]->Write(); | |
411 | } | |
412 | ||
413 | */ | |
414 | ||
415 | if(simflag==false){ | |
416 | fEventStat->Write(); | |
417 | fEvent->Write(); | |
418 | fdEdXvsPt->Write(); | |
419 | fdEdXnSigmaElecvsPt->Write(); | |
420 | fTOFbetavsPt->Write(); | |
421 | fTOFnSigmaElecvsPt->Write(); | |
422 | fdEdXvsPtWithCut->Write(); | |
423 | fPtElec[0]->Write(); | |
424 | fPtElec[1]->Write(); | |
425 | } | |
426 | fNelc_pos->Write(); | |
427 | fNelc_all->Write(); | |
428 | for(int i=0;i<10;i++){ | |
429 | fNelc_pos_cent[i]->Write(); | |
430 | fNelc_all_cent[i]->Write(); | |
431 | } | |
432 | fNelc_all_pT->Write(); | |
433 | //d_tree->Write(); | |
434 | hCentrality->Write(); | |
435 | hV0AC->Write(); | |
436 | hV0AC_Ntrk->Write(); | |
437 | hV0AC_NaccTrcklts->Write(); | |
438 | for(int j=0;j<11;j++){ | |
439 | for(int i=0;i<7;i++){ | |
440 | hmasspt[i][j]->Write(); | |
441 | hmasspt_weight[i][j]->Write(); | |
442 | } | |
443 | } | |
444 | ||
445 | //d_ntpair->Write(); | |
446 | fout->Close(); | |
447 | } | |
448 | ||
449 | //____________________________________________ | |
450 | void ana_sgl::loop_a_file(char *file){ | |
451 | ||
452 | TFile *treefile = TFile::Open(file); | |
453 | TDirectory *d = (TDirectory*)treefile->Get("PWG3_dielectron"); | |
454 | if(d==0){ | |
455 | cout<<" PWG3_dielectron is not found "<<endl; | |
456 | return ; | |
457 | } | |
458 | ||
459 | //TTree *tree = (TTree*)d->Get("t"); | |
460 | TTree *tree = (TTree*)d->Get("tree_MultiDie_CENT1"); | |
461 | if(tree == 0) { | |
462 | cout <<"tree is not found in "<<file<<endl; | |
463 | treefile->Close(); | |
464 | return; | |
465 | } | |
466 | cout << file <<" is opened"<<endl; | |
467 | if(simflag==false){ | |
468 | //add_histograms(treefile); | |
469 | } | |
470 | Init(tree); | |
471 | Loop(); | |
472 | tree->Clear(); | |
473 | d->Clear(); | |
474 | tree->Delete(); | |
475 | d->Delete(); | |
476 | delete d; | |
477 | treefile->Close(); | |
478 | delete treefile; | |
479 | ||
480 | ||
481 | ||
482 | cout <<"one file processed"<<endl; | |
483 | } | |
484 | ||
485 | //____________________________________________ | |
486 | void ana_sgl::ana_event(int ientry, int jentry){ | |
487 | ||
488 | //select trigger class: | |
489 | if(sel_trigger==2){ | |
490 | if(kTriggerCent<100){ | |
491 | return ; | |
492 | } | |
493 | }else if(sel_trigger==1){ | |
494 | if(kTriggerCent%100<10){ | |
495 | return ; | |
496 | } | |
497 | }else if(sel_trigger==0){ | |
498 | if(kTriggerCent%10!=1){ | |
499 | return ; | |
500 | } | |
501 | } | |
502 | ||
503 | // if(fkRunNumber>138350){ | |
504 | /* | |
505 | if(kMag>0){ | |
506 | magnetic_field_mm = false; | |
507 | }else{ | |
508 | magnetic_field_mm = true; | |
509 | } | |
510 | */ | |
511 | if(fkRunNumber>169591){ | |
512 | magnetic_field_mm = false; | |
513 | }else{ | |
514 | magnetic_field_mm = true; | |
515 | } | |
516 | ||
517 | if(ientry%1000==0){ | |
518 | cout<<" event processing "<<ientry<<" / "<<d_evt<<" / "<<d_event<<" : trigger "<<kTriggerCent<<endl; | |
519 | } | |
520 | ||
521 | d_nelc = 0; | |
522 | for(int i=0;i<10;i++){ | |
523 | nelec_pos[i] = 0; | |
524 | } | |
525 | ||
526 | hCentrality->Fill(fkCentrality); | |
527 | hV0AC->Fill(fkV0C, fkV0A); | |
528 | hV0AC_Ntrk->Fill(fkNTrk, fkV0C+fkV0A); | |
529 | hV0AC_NaccTrcklts->Fill(fkNaccTrcklts, fkV0C+fkV0A); | |
530 | fill_to_tree_variables(); | |
531 | ||
532 | for(int i=0;i<fkNPar;i++){ | |
533 | ||
534 | if(simflag==false){ | |
535 | ||
536 | //fdEdXvsPtWithCut->Fill(kPIn[i], kTPCsignal[i]); | |
537 | ||
538 | if(GlobalTrackcut(i)==false){ | |
539 | continue; | |
540 | } | |
541 | fill_histograms(i); | |
542 | } | |
543 | ||
544 | ||
545 | fill_to_tree_track_variables(i); | |
546 | ||
547 | //if(PairTrackcut(i)==false){ | |
548 | //continue; | |
549 | //} | |
550 | ||
551 | if(fkCentrality>0 && fkCentrality<10){ | |
552 | int iptbin = (int)(kP[i]/0.5); | |
553 | if(iptbin>=10){ | |
554 | iptbin=9; | |
555 | } | |
556 | nelec_pos[iptbin]++; | |
557 | } | |
558 | ||
559 | if(kCharge[i]>0){ | |
560 | etrk e = etrk( | |
561 | fkCentrality, fkXvPrim, fkYvPrim, fkZvPrim, | |
562 | kXv[i], kYv[i], kZv[i], | |
563 | kPx[i], kPy[i], kPz[i], kPt[i], | |
564 | kEta[i], kPhi[i], kTheta[i], | |
565 | kTPCsignal[i], kTOFbeta[i], | |
566 | kE[i], kDeltaPhi[i], kDeltaEta[i], | |
567 | kTPCnSigmaEle[i], kTPCnSigmaPio[i], kTPCnSigmaKao[i], kTPCnSigmaPro[i], | |
568 | kTOFnSigmaEle[i], kTOFnSigmaPio[i], kTOFnSigmaKao[i], kTOFnSigmaPro[i], | |
569 | kITSsignal[i], kNclsITS[i], kNclsTPC[i], kLegDistXY[i], kLegDist[i] | |
570 | ); | |
571 | vep_tmp.push_back(e); | |
572 | }else{ | |
573 | etrk e = etrk( | |
574 | fkCentrality, fkXvPrim, fkYvPrim, fkZvPrim, | |
575 | kXv[i], kYv[i], kZv[i], | |
576 | kPx[i], kPy[i], kPz[i], kPt[i], | |
577 | kEta[i], kPhi[i], kTheta[i], | |
578 | kTPCsignal[i], kTOFbeta[i], | |
579 | kE[i], kDeltaPhi[i], kDeltaEta[i], | |
580 | kTPCnSigmaEle[i], kTPCnSigmaPio[i], kTPCnSigmaKao[i], kTPCnSigmaPro[i], | |
581 | kTOFnSigmaEle[i], kTOFnSigmaPio[i], kTOFnSigmaKao[i], kTOFnSigmaPro[i], | |
582 | kITSsignal[i], kNclsITS[i], kNclsTPC[i], kLegDistXY[i], kLegDist[i] | |
583 | ); | |
584 | vem_tmp.push_back(e); | |
585 | } | |
586 | } | |
587 | ||
588 | //////// fill to the tree ////////////////// | |
589 | //d_tree->Fill(); | |
590 | ||
591 | ||
592 | ||
593 | if(d_conv_flag==true){ | |
594 | check_conversion_pairs(vep_tmp, vem_tmp); | |
595 | } | |
596 | ||
597 | check_ghost_pairs(vep_tmp); | |
598 | check_ghost_pairs(vem_tmp); | |
599 | randomize_pool(vep_tmp, vem_tmp); | |
600 | if(d_conv_flag==false){ | |
601 | check_conversion_pairs(vep, vem); | |
602 | } | |
603 | ||
604 | fNelc_pos->Fill(vep.size(), vem.size()); | |
605 | fNelc_all->Fill(vep.size()+vem.size()); | |
606 | ||
607 | if(fkCentrality>0 && fkCentrality<10){ | |
608 | for(int i=0;i<10;i++){ | |
609 | fNelc_all_pT->Fill(0.25+0.5*i, nelec_pos[i]); | |
610 | } | |
611 | } | |
612 | ||
613 | ||
614 | calc_pair(vep, vem); | |
615 | ||
616 | int icent = (int)(fkCentrality/10.0); | |
617 | int izbin = (int)((fkZvPrim+10)/2.0); | |
618 | if(icent<0) icent=0; | |
619 | if(icent>=NCENT) icent=NCENT-1; | |
620 | if(izbin<0) izbin=0; | |
621 | if(izbin>=NZBIN) izbin=NZBIN-1; | |
622 | ||
623 | fNelc_pos_cent[icent]->Fill(vep.size(), vem.size()); | |
624 | fNelc_all_cent[icent]->Fill(vep.size()+vem.size()); | |
625 | ||
626 | ||
627 | ||
628 | d_vep[d_ibuf[izbin][icent]][izbin][icent].clear(); | |
629 | vector<etrk>::iterator iep; | |
630 | for(iep = vep.begin();iep != vep.end();++iep) { | |
631 | d_vep[d_ibuf[izbin][icent]][izbin][icent].push_back(*iep); | |
632 | d_poolp[izbin][icent].push_back(*iep); | |
633 | if(d_poolp[izbin][icent].size()>MAXPOOL) { | |
634 | d_poolp[izbin][icent].pop_front(); | |
635 | } | |
636 | } | |
637 | d_vem[d_ibuf[izbin][icent]][izbin][icent].clear(); | |
638 | vector<etrk>::iterator iem; | |
639 | for(iem = vem.begin();iem != vem.end();++iem) { | |
640 | d_vem[d_ibuf[izbin][icent]][izbin][icent].push_back(*iem); | |
641 | d_poolm[izbin][icent].push_back(*iem); | |
642 | if(d_poolm[izbin][icent].size()>MAXPOOL) { | |
643 | d_poolm[izbin][icent].pop_front(); | |
644 | } | |
645 | } | |
646 | // Update the buffer pointer | |
647 | d_ibuf[izbin][icent]++; | |
648 | if(d_ibuf[izbin][icent]>= NBUF) d_ibuf[izbin][icent]=0; | |
649 | ||
650 | /////////////////////////////////////////// | |
651 | vem.clear(); | |
652 | vep.clear(); | |
653 | vem_tmp.clear(); | |
654 | vep_tmp.clear(); | |
655 | d_evt ++; | |
656 | d_event ++; | |
657 | ||
658 | } | |
659 | ||
660 | ||
661 | //____________________________________________ | |
662 | bool ana_sgl::kTOFcut(int itrk){ | |
663 | ||
664 | if(kTOFsignal[itrk]!=9999 && kTOFbeta[itrk]>0.3 && | |
665 | (kTOFnSigmaEle[itrk]<3 && kTOFnSigmaEle[itrk]>-3) | |
666 | /* | |
667 | (kTOFnSigmaKao>3 || kTOFnSigmaKao<-3) && | |
668 | (kTOFnSigmaPro>3 || kTOFnSigmaPro<-3) | |
669 | */ | |
670 | ||
671 | ){ | |
672 | ||
673 | return true; | |
674 | }else{ | |
675 | return false; | |
676 | } | |
677 | } | |
678 | ||
679 | //____________________________________________ | |
680 | bool ana_sgl::GlobalTrackcut(int itrk){ | |
681 | ||
682 | ||
683 | if(( | |
684 | (d_flag_kaon_veto == true && (kTPCnSigmaKao[itrk]<d_dedx_kaon_veto_low || kTPCnSigmaKao[itrk]>d_dedx_kaon_veto_high)) || | |
685 | (d_flag_kaon_veto == false) | |
686 | ) | |
687 | && ( | |
688 | (d_flag_proton_veto == true && (kTPCnSigmaPro[itrk]<d_dedx_proton_veto_low || kTPCnSigmaPro[itrk]>d_dedx_proton_veto_high)) || | |
689 | (d_flag_proton_veto == false) | |
690 | ) | |
691 | ){ | |
692 | fdEdXvsPtWithCut->Fill(kPIn[itrk], kTPCsignal[itrk]); | |
693 | } | |
694 | ||
695 | if( kNclsTPC[itrk]>120 | |
696 | && kTPCsignal[itrk]>d_tpc_dedx_low && kTPCsignal[itrk]<d_tpc_dedx_high | |
697 | && ( | |
698 | (d_flag_tof_cut==true && kTOFnSigmaEle[itrk]>d_tof_low && kTOFnSigmaEle[itrk]<d_tof_high) || | |
699 | (d_flag_tof_cut==false) | |
700 | ) | |
701 | && ( | |
702 | (d_flag_pt_cut == true && kPt[itrk]>d_pt_cut_low && kPt[itrk]<d_pt_cut_high) || | |
703 | (d_flag_pt_cut == false) | |
704 | ) | |
705 | && ( | |
706 | (d_flag_kaon_veto == true && (kTPCnSigmaKao[itrk]<d_dedx_kaon_veto_low || kTPCnSigmaKao[itrk]>d_dedx_kaon_veto_high)) || | |
707 | (d_flag_kaon_veto == false) | |
708 | ) | |
709 | && ( | |
710 | (d_flag_proton_veto == true && (kTPCnSigmaPro[itrk]<d_dedx_proton_veto_low || kTPCnSigmaPro[itrk]>d_dedx_proton_veto_high)) || | |
711 | (d_flag_proton_veto == false) | |
712 | ) | |
713 | ){ | |
714 | return true; | |
715 | }else{ | |
716 | return false; | |
717 | } | |
718 | } | |
719 | ||
720 | /* | |
721 | //____________________________________________ | |
722 | bool ana_sgl::PairTrackcut(int itrk){ | |
723 | ||
724 | if( (kTOFnSigmaEle[itrk]<3 && kTOFnSigmaEle[itrk]>-3) && | |
725 | (kTOFnSigmaKao[itrk]>3||kTOFnSigmaKao[itrk]<-3) && | |
726 | (kTOFnSigmaPro[itrk]>3||kTOFnSigmaPro[itrk]<-3) //&& | |
727 | //(kTPCnSigmaEle[itrk]-1.65226*exp(-kPIn[itrk]*kPIn[itrk]*1.60890)+0.838)<1 && | |
728 | //(kTPCnSigmaEle[itrk]-1.65226*exp(-kPIn[itrk]*kPIn[itrk]*1.60890)+0.838)>-1 | |
729 | ){ | |
730 | return true; | |
731 | }else{ | |
732 | return false; | |
733 | } | |
734 | ||
735 | ||
736 | // if( | |
737 | // !(kTPCsignal>75-210*(kPIn-0.5) && | |
738 | // kTPCsignal<100-190*(kPIn-0.5)) && | |
739 | // !(kTPCsignal>130-120*(kPIn-0.5) && | |
740 | // kTPCsignal<135-80*(kPIn-0.5)) && | |
741 | // (kTPCnSigmaEle-1.65226*exp(-kPIn*kPIn*1.60890)+0.838)<1 && | |
742 | // (kTPCnSigmaEle-1.65226*exp(-kPIn*kPIn*1.60890)+0.838)>-1 | |
743 | // ){ | |
744 | // return true; | |
745 | // }else{ | |
746 | // return false; | |
747 | // } | |
748 | ||
749 | ||
750 | } | |
751 | */ | |
752 | ||
753 | //____________________________________________ | |
754 | void ana_sgl::fill_histograms(int itrk){ | |
755 | /* | |
756 | hdedx_pt->Fill(kPIn[itrk], kTPCsignal[itrk]); | |
757 | hbetatof_pt->Fill(kPIn[itrk], kTOFbeta[itrk]); | |
758 | ||
759 | if(kTOFcut(itrk)==true){ | |
760 | hdedx_tof_pt->Fill(kPIn[itrk], kTPCsignal[itrk]); | |
761 | hbetatof_tof_pt->Fill(kPIn[itrk], kTOFbeta[itrk]); | |
762 | } | |
763 | ||
764 | int iptbin = (int)((kPIn[itrk])/fBinWidth); | |
765 | if(iptbin>=0 && iptbin<nHistos){ | |
766 | hdedx[iptbin]->Fill(kTPCsignal[itrk]); | |
767 | if(kTOFcut(itrk)==true)hdedx_tof[iptbin]->Fill(kTPCsignal[itrk]); | |
768 | } | |
769 | */ | |
770 | ||
771 | if(kCharge[itrk]>0){ | |
772 | fPtElec[1]->Fill(kPIn[itrk]); | |
773 | }else{ | |
774 | fPtElec[0]->Fill(kPIn[itrk]); | |
775 | } | |
776 | ||
777 | hdedx_pt->Fill(kPIn[itrk], kTPCsignal[itrk]); | |
778 | hbetatof_pt->Fill(kPIn[itrk], kTOFbeta[itrk]); | |
779 | ||
780 | if(kTOFnSigmaEle[itrk]>-3 && kTOFnSigmaEle[itrk]<3){ | |
781 | hdedx_tof_elec_pt->Fill(kPIn[itrk], kTPCsignal[itrk]); | |
782 | hbetatof_tof_elec_pt->Fill(kPIn[itrk], kTOFbeta[itrk]); | |
783 | ||
784 | if(kE[itrk]/kPIn[itrk]>0.7 && kE[itrk]/kPIn[itrk]<1.3){ | |
785 | hdedx_tof_elec_emc_pt->Fill(kPIn[itrk], kTPCsignal[itrk]); | |
786 | hbetatof_emc_pt->Fill(kPIn[itrk], kTOFbeta[itrk]); | |
787 | } | |
788 | ||
789 | if( (kTOFnSigmaKao[itrk]>3 || kTOFnSigmaKao[itrk]<-3) && | |
790 | (kTOFnSigmaPro[itrk]>3 || kTOFnSigmaPro[itrk]<-3) | |
791 | ){ | |
792 | hdedx_tof_all_pt->Fill(kPIn[itrk], kTPCsignal[itrk]); | |
793 | hbetatof_tof_all_pt->Fill(kPIn[itrk], kTPCsignal[itrk]); | |
794 | if(kE[itrk]/kPIn[itrk]>0.7 && kE[itrk]/kPIn[itrk]<1.3){ | |
795 | hdedx_tof_all_emc_pt->Fill(kPIn[itrk], kTPCsignal[itrk]); | |
796 | hbetatof_tof_all_emc_pt->Fill(kPIn[itrk], kTPCsignal[itrk]); | |
797 | } | |
798 | } | |
799 | } | |
800 | ||
801 | if(kE[itrk]/kPIn[itrk]>0.7 && kE[itrk]/kPIn[itrk]<1.3){ | |
802 | hdedx_emc_pt->Fill(kPIn[itrk], kTPCsignal[itrk]); | |
803 | hbetatof_emc_pt->Fill(kPIn[itrk], kTPCsignal[itrk]); | |
804 | } | |
805 | ||
806 | ||
807 | ||
808 | ||
809 | int iptbin = (int)((kPIn[itrk])/fBinWidth); | |
810 | if(iptbin>=0 && iptbin<nHistos){ | |
811 | hdedx[iptbin]->Fill(kTPCsignal[itrk]); | |
812 | //if(kTOFcut(itrk)==true)hdedx_tof[iptbin]->Fill(kTPCsignal[itrk]); | |
813 | ||
814 | if(kTOFnSigmaEle[itrk]>-3 && kTOFnSigmaEle[itrk]<3){ | |
815 | hdedx_tof_elec[iptbin]->Fill(kTPCsignal[itrk]); | |
816 | ||
817 | if(kE[itrk]/kPIn[itrk]>0.7 && kE[itrk]/kPIn[itrk]<1.3){ | |
818 | hdedx_tof_elec_emc[iptbin]->Fill(kTPCsignal[itrk]); | |
819 | } | |
820 | if( (kTOFnSigmaKao[itrk]>3 || kTOFnSigmaKao[itrk]<-3) && | |
821 | (kTOFnSigmaPro[itrk]>3 || kTOFnSigmaPro[itrk]<-3) | |
822 | ){ | |
823 | hdedx_tof_all[iptbin]->Fill(kTPCsignal[itrk]); | |
824 | if(kE[itrk]/kPIn[itrk]>0.7 && kE[itrk]/kPIn[itrk]<1.3){ | |
825 | hdedx_tof_all_emc[iptbin]->Fill(kTPCsignal[itrk]); | |
826 | } | |
827 | } | |
828 | } | |
829 | if(kE[itrk]/kPIn[itrk]>0.7 && kE[itrk]/kPIn[itrk]<1.3){ | |
830 | hdedx_emc[iptbin]->Fill(kTPCsignal[itrk]); | |
831 | } | |
832 | } | |
833 | ||
834 | } | |
835 | ||
836 | //____________________________________________ | |
837 | void ana_sgl::add_histograms(TFile *fin){ | |
838 | ||
839 | TDirectory *d = (TDirectory*)fin->Get("PWG3_dielectron"); | |
840 | if(d==0){ | |
841 | cout<<" PWG3_dielectron is not found "<<endl; | |
842 | return ; | |
843 | } | |
844 | fEventStat->Add((TH1D*)d->Get("hEventStat_MultiDie_CENT1")); | |
845 | TList *list = (TList*)d->Get("jpsi_QA_CENT1"); | |
846 | TList *listQA = (TList*)list->FindObject("QAElectron"); | |
847 | fEvent->Add((TH1D*)listQA->FindObject("Event")); | |
848 | fdEdXvsPt->Add((TH1D*)listQA->FindObject("dEdXvsPt")); | |
849 | fdEdXnSigmaElecvsPt->Add((TH1D*)listQA->FindObject("fdEdXnSigmaElecvsPt")); | |
850 | fTOFbetavsPt->Add((TH1D*)listQA->FindObject("fTOFbetavsPt")); | |
851 | fTOFnSigmaElecvsPt->Add((TH1D*)listQA->FindObject("fTOFnSigmaElecvsPt")); | |
852 | ||
853 | cout<<" ana_sgl::add_histograms "<<fin->GetName()<<" done "<<endl; | |
854 | ||
855 | } | |
856 | ||
857 | //_____________________________________________ | |
858 | void ana_sgl::fill_to_tree_variables(void){ | |
859 | ||
860 | d_cent = fkCentrality; | |
861 | d_ntrk= fkNTrk; | |
862 | d_xvprim= fkXvPrim; | |
863 | d_yvprim= fkYvPrim; | |
864 | d_zvprim= fkZvPrim; | |
865 | d_nacctrklets= fkNaccTrcklts; | |
866 | d_xres= fkXRes; | |
867 | d_yres= fkYRes; | |
868 | d_zres= fkZRes; | |
869 | ||
870 | } | |
871 | //_____________________________________________ | |
872 | void ana_sgl::fill_to_tree_track_variables(int itrk){ | |
873 | ||
874 | ||
875 | d_px[d_nelc]= kPx[itrk]; | |
876 | d_py[d_nelc]= kPy[itrk]; | |
877 | d_pz[d_nelc]= kPz[itrk]; | |
878 | d_p[d_nelc]= kP[itrk]; | |
879 | d_pt[d_nelc]= kPt[itrk]; | |
880 | d_xv[d_nelc]= kXv[itrk]; | |
881 | d_yv[d_nelc]= kYv[itrk]; | |
882 | d_zv[d_nelc]= kZv[itrk]; | |
883 | d_phi[d_nelc]= kPhi[itrk]; | |
884 | d_theta[d_nelc]= kTheta[itrk]; | |
885 | d_eta[d_nelc]= kEta[itrk]; | |
886 | d_c[d_nelc]= kCharge[itrk]; | |
887 | d_nclusITS[d_nelc]= kNclsITS[itrk]; | |
888 | d_nclusTPC[d_nelc]= kNclsTPC[itrk]; | |
889 | d_nclusTPCiter[d_nelc]= kNclsTPCiter1[itrk]; | |
890 | d_nfclusTPC[d_nelc]= kNFclsTPC[itrk]; | |
891 | d_nfclusTPCr[d_nelc]= kNFclsTPCr[itrk]; | |
892 | d_nfclusTPCrFrac[d_nelc]= kNFclsTPCrFrac[itrk]; | |
893 | d_TPCsignalN[d_nelc]= kTPCsignalN[itrk]; | |
894 | d_TPCsignalNfrac[d_nelc]= kTPCsignalNfrac[itrk]; | |
895 | d_TPCchi2cl[d_nelc]= kTPCchi2Cl[itrk]; | |
896 | d_trkstat[d_nelc]= kTrackStatus[itrk]; | |
897 | d_nclsTRD[d_nelc]= kNclsTRD[itrk]; | |
898 | d_TRDntracklets[d_nelc]= kTRDntracklets[itrk]; | |
899 | d_TRDpidquality[d_nelc]= kTRDpidQuality[itrk]; | |
900 | d_TRDprobEle[d_nelc]= kTRDprobEle[itrk]; | |
901 | d_TRDprobPio[d_nelc]= kTRDprobPio[itrk]; | |
902 | d_impactXY[d_nelc]= kImpactParXY[itrk]; | |
903 | d_impactZ[d_nelc]= kImpactParZ[itrk]; | |
904 | d_tracklength[d_nelc]= kTrackLength[itrk]; | |
905 | d_ITSsignal[d_nelc]= kITSsignal[itrk]; | |
906 | d_ITSnsigmaEle[d_nelc]= kITSnSigmaEle[itrk]; | |
907 | d_ITSnsigmaPio[d_nelc]= kITSnSigmaPio[itrk]; | |
908 | d_ITSnsigmaMuo[d_nelc]= kITSnSigmaMuo[itrk]; | |
909 | d_ITSnsigmaKao[d_nelc]= kITSnSigmaKao[itrk]; | |
910 | d_ITSnsigmaPro[d_nelc]= kITSnSigmaPro[itrk]; | |
911 | d_PIn[d_nelc]= kPIn[itrk]; | |
912 | d_TPCsignal[d_nelc]= kTPCsignal[itrk]; | |
913 | d_TOFsignal[d_nelc]= kTOFsignal[itrk]; | |
914 | d_TOFbeta[d_nelc]= kTOFbeta[itrk]; | |
915 | d_TPCnSigmaEle[d_nelc]= kTPCnSigmaEle[itrk]; | |
916 | d_TPCnSigmaPio[d_nelc]= kTPCnSigmaPio[itrk]; | |
917 | d_TPCnSigmaMuo[d_nelc]= kTPCnSigmaMuo[itrk]; | |
918 | d_TPCnSigmaKao[d_nelc]= kTPCnSigmaKao[itrk]; | |
919 | d_TPCnSigmaPro[d_nelc]= kTPCnSigmaPro[itrk]; | |
920 | d_TOFnSigmaEle[d_nelc]= kTOFnSigmaEle[itrk]; | |
921 | d_TOFnSigmaPio[d_nelc]= kTOFnSigmaPio[itrk]; | |
922 | d_TOFnSigmaMuo[d_nelc]= kTOFnSigmaMuo[itrk]; | |
923 | d_TOFnSigmaKao[d_nelc]= kTOFnSigmaKao[itrk]; | |
924 | d_TOFnSigmaPro[d_nelc]= kTOFnSigmaPro[itrk]; | |
925 | ||
926 | d_chi2ndf[d_nelc]= kChi2NDF[itrk]; | |
927 | d_E[d_nelc] = kE[itrk]; | |
928 | d_dphi[d_nelc] = kDeltaPhi[itrk]; | |
929 | d_deta[d_nelc] = kDeltaEta[itrk]; | |
930 | ||
931 | d_nelc++; | |
932 | ||
933 | } | |
934 | ||
935 | //____________________________________________________________________ | |
936 | void ana_sgl::randomize_pool(vector<etrk> e1, vector<etrk> e2){ | |
937 | ||
938 | int size1 = e1.size(); | |
939 | int used_index[1000]; | |
940 | for(int i=0;i<1000;i++){ | |
941 | used_index[i] = -1; | |
942 | } | |
943 | for(int i=0;i<size1;i++){ | |
944 | used_index[i] = 0; | |
945 | } | |
946 | ||
947 | for(unsigned int i=0;i<size1;i++){ | |
948 | int j = (int)(gRandom->Uniform(0,size1)); | |
949 | while(used_index[j]==1){ | |
950 | j = (int)(gRandom->Uniform(0,size1)); | |
951 | } | |
952 | if( (e1[j].ghost_flag==1) && | |
953 | ( | |
954 | (d_conv_flag==true && e1[j].conv_flag==1) || | |
955 | (d_conv_flag==false) | |
956 | ) | |
957 | ){ | |
958 | vep.push_back(e1[j]); | |
959 | } | |
960 | used_index[j] = 1; | |
961 | } | |
962 | ||
963 | ||
964 | int size2 = e2.size(); | |
965 | for(int i=0;i<1000;i++){ | |
966 | used_index[i] = -1; | |
967 | } | |
968 | for(int i=0;i<size2;i++){ | |
969 | used_index[i] = 0; | |
970 | } | |
971 | ||
972 | for(unsigned int i=0;i<size2;i++){ | |
973 | int j = (int)(gRandom->Uniform(0,size2)); | |
974 | while(used_index[j]==1){ | |
975 | j = (int)(gRandom->Uniform(0,size2)); | |
976 | } | |
977 | if( (e2[j].ghost_flag==1) && | |
978 | ( | |
979 | (d_conv_flag==true && e2[j].conv_flag==1) || | |
980 | (d_conv_flag==false) | |
981 | ) | |
982 | ){ | |
983 | vem.push_back(e2[j]); | |
984 | } | |
985 | used_index[j] = 1; | |
986 | } | |
987 | } | |
988 | ||
989 | //____________________________________________________________________ | |
990 | void ana_sgl::calc_pair(vector<etrk> vep, vector<etrk> vem){ | |
991 | ||
992 | //vector<etrk>::iterator iep, iem; | |
993 | ||
994 | vector<etrk>::iterator iep; | |
995 | vector<etrk>::iterator iem; | |
996 | ||
997 | //cout<<vep.size()<<" "<<vem.size()<<endl; | |
998 | ||
999 | ///unlike pairs | |
1000 | for(iep=vep.begin(); iep!=vep.end(); ++iep){ | |
1001 | for(iem=vem.begin(); iem!=vem.end(); ++iem){ | |
1002 | if(PairTrackcut(iep, iem)==false) continue; | |
1003 | fill_pair(iep, iem, 0); | |
1004 | } | |
1005 | } | |
1006 | for(iep=vep.begin(); iep!=vep.end(); ++iep){ | |
1007 | vector<etrk>::iterator iep2=iep; | |
1008 | ++iep2; | |
1009 | for(iem=iep2; iem!=vep.end(); ++iem){ | |
1010 | if(PairTrackcut(iep, iem)==false) continue; | |
1011 | fill_pair(iep, iem, 1); | |
1012 | } | |
1013 | } | |
1014 | ||
1015 | for(iep=vem.begin(); iep!=vem.end(); ++iep){ | |
1016 | vector<etrk>::iterator iep2=iep; | |
1017 | ++iep2; | |
1018 | for(iem=iep2; iem!=vem.end(); ++iem){ | |
1019 | if(PairTrackcut(iep, iem)==false) continue; | |
1020 | fill_pair(iep, iem, 2); | |
1021 | } | |
1022 | } | |
1023 | ||
1024 | int icent = (int)(fkCentrality/10.0); | |
1025 | int izbin = (int)((fkZvPrim+10)/2.0); | |
1026 | if(icent<0) icent=0; | |
1027 | if(icent>=NCENT) icent=NCENT-1; | |
1028 | if(izbin<0) izbin=0; | |
1029 | if(izbin>=NZBIN) izbin=NZBIN-1; | |
1030 | ||
1031 | ||
1032 | int nmixed; | |
1033 | if(vep.size()>0) { | |
1034 | // | |
1035 | // Now mixed event for +- pairs | |
1036 | // | |
1037 | nmixed = 0; | |
1038 | for(int ibuf=0;(nmixed<NMix);ibuf++) { | |
1039 | int ntry = 0; | |
1040 | while(ntry<MAX_TRY) { | |
1041 | reshuffle_buffer(d_vem[ibuf][izbin][icent],d_poolm[izbin][icent]); | |
1042 | ntry++; | |
1043 | } | |
1044 | for(iep=vep.begin(); iep!=vep.end(); ++iep){ | |
1045 | for(iem=d_vem[ibuf][izbin][icent].begin(); | |
1046 | iem!=d_vem[ibuf][izbin][icent].end(); ++iem){ | |
1047 | if(PairTrackcut(iep, iem)==false) continue; | |
1048 | fill_pair(iep,iem,3); | |
1049 | } | |
1050 | } | |
1051 | ++nmixed; | |
1052 | }//for(ibuf) | |
1053 | } | |
1054 | ||
1055 | if(vem.size()>0) { | |
1056 | // | |
1057 | // Now mixed event for +- pairs | |
1058 | // | |
1059 | nmixed = 0; | |
1060 | for(int ibuf=0;(nmixed<NMix);ibuf++) { | |
1061 | int ntry = 0; | |
1062 | while(ntry<MAX_TRY) { | |
1063 | reshuffle_buffer(d_vep[ibuf][izbin][icent],d_poolp[izbin][icent]); | |
1064 | ntry++; | |
1065 | } | |
1066 | for(iem=vem.begin(); iem!=vem.end(); ++iem){ | |
1067 | for(iep=d_vep[ibuf][izbin][icent].begin(); | |
1068 | iep!=d_vep[ibuf][izbin][icent].end(); ++iep){ | |
1069 | if(PairTrackcut(iep, iem)==false) continue; | |
1070 | fill_pair(iep,iem,4); | |
1071 | } | |
1072 | } | |
1073 | ++nmixed; | |
1074 | }//for(ibuf) | |
1075 | } | |
1076 | ||
1077 | ||
1078 | if(vep.size()>0) { | |
1079 | // | |
1080 | // Now mixed event for ++ pairs | |
1081 | // | |
1082 | nmixed = 0; | |
1083 | for(int ibuf=0;(nmixed<NMix);ibuf++) { | |
1084 | int ntry = 0; | |
1085 | while(ntry<MAX_TRY) { | |
1086 | reshuffle_buffer(d_vep[ibuf][izbin][icent],d_poolp[izbin][icent]); | |
1087 | ntry++; | |
1088 | } | |
1089 | for(iep=vep.begin(); iep!=vep.end(); ++iep){ | |
1090 | for(iem=d_vep[ibuf][izbin][icent].begin(); | |
1091 | iem!=d_vep[ibuf][izbin][icent].end(); ++iem){ | |
1092 | if(PairTrackcut(iep, iem)==false) continue; | |
1093 | fill_pair(iep,iem,5); | |
1094 | } | |
1095 | } | |
1096 | ++nmixed; | |
1097 | }//for(ibuf) | |
1098 | } | |
1099 | ||
1100 | if(vem.size()>0) { | |
1101 | // | |
1102 | // Now mixed event for +- pairs | |
1103 | // | |
1104 | nmixed = 0; | |
1105 | for(int ibuf=0;(nmixed<NMix);ibuf++) { | |
1106 | int ntry = 0; | |
1107 | while(ntry<MAX_TRY) { | |
1108 | reshuffle_buffer(d_vem[ibuf][izbin][icent],d_poolm[izbin][icent]); | |
1109 | ntry++; | |
1110 | } | |
1111 | for(iem=vem.begin(); iem!=vem.end(); ++iem){ | |
1112 | for(iep=d_vem[ibuf][izbin][icent].begin(); | |
1113 | iep!=d_vem[ibuf][izbin][icent].end(); ++iep){ | |
1114 | if(PairTrackcut(iep, iem)==false) continue; | |
1115 | fill_pair(iep,iem,6); | |
1116 | } | |
1117 | } | |
1118 | ++nmixed; | |
1119 | }//for(ibuf) | |
1120 | } | |
1121 | } | |
1122 | ||
1123 | ||
1124 | //___________________________________________________________________________ | |
1125 | //void ana_sgl::fill_pair(etrk* iep, etrk* iem, int type){ | |
1126 | //void ana_sgl::fill_pair(etrk iep, etrk iem, int type){ | |
1127 | void ana_sgl::fill_pair(vector<etrk>::iterator iep, vector<etrk>::iterator iem, int type){ | |
1128 | ||
1129 | ||
1130 | d_pairtype = type; | |
1131 | ||
1132 | calc_vars(iep, iem, d_mass, d_phiv, d_pxpair, d_pypair, d_pzpair, | |
1133 | d_ptpair, d_epair, d_phipair, d_etapair, d_cos, d_psi); | |
1134 | ||
1135 | if(type==0||type==1||type==2||type==3||type==5){ | |
1136 | d_centrality = iep->cent; | |
1137 | d_prim_xv = iep->pxv; | |
1138 | d_prim_yv = iep->pyv; | |
1139 | d_prim_zv = iep->pzv; | |
1140 | }else if(type==4 || type==6){ | |
1141 | d_centrality = iem->cent; | |
1142 | d_prim_xv = iem->pxv; | |
1143 | d_prim_yv = iem->pyv; | |
1144 | d_prim_zv = iem->pzv; | |
1145 | ||
1146 | } | |
1147 | ||
1148 | ||
1149 | //cout<<iep->px<<" "<<iem->px<<" "<<d_centrality<<endl; | |
1150 | ||
1151 | d_cent1 = iep->cent; | |
1152 | d_xv1 = iep->xv; | |
1153 | d_yv1 = iep->yv; | |
1154 | d_zv1 = iep->zv; | |
1155 | d_px1 = iep->px; | |
1156 | d_py1 = iep->py; | |
1157 | d_pz1 = iep->pz; | |
1158 | d_pt1 = iep->pt; | |
1159 | d_eta1 = iep->eta; | |
1160 | d_phi1 = iep->phi; | |
1161 | d_theta1 = iep->theta; | |
1162 | d_tpc1 = iep->tpc; | |
1163 | d_ntpc_ele1 = iep->ntpc_ele; | |
1164 | d_ntpc_pio1 = iep->ntpc_pio; | |
1165 | d_ntpc_kao1 = iep->ntpc_kao; | |
1166 | d_ntpc_pro1 = iep->ntpc_pro; | |
1167 | d_ntof_ele1 = iep->ntof_ele; | |
1168 | d_ntof_pio1 = iep->ntof_pio; | |
1169 | d_ntof_kao1 = iep->ntof_kao; | |
1170 | d_ntof_pro1 = iep->ntof_pro; | |
1171 | d_its1 = iep->its; | |
1172 | d_nits1 = iep->nits; | |
1173 | d_ntpc1 = iep->ntpc; | |
1174 | d_e1 = iep->e; | |
1175 | d_dphi1 = iep->dphi; | |
1176 | d_deta1 = iep->deta; | |
1177 | d_dcaxy1 = iep->dcaxy; | |
1178 | d_dcaz1 = iep->dcaz; | |
1179 | d_conv1 = iep->conv_flag; | |
1180 | ||
1181 | d_cent2 = iem->cent; | |
1182 | d_xv2 = iem->xv; | |
1183 | d_yv2 = iem->yv; | |
1184 | d_zv2 = iem->zv; | |
1185 | d_px2 = iem->px; | |
1186 | d_py2 = iem->py; | |
1187 | d_pz2 = iem->pz; | |
1188 | d_pt2 = iem->pt; | |
1189 | d_eta2 = iem->eta; | |
1190 | d_phi2 = iem->phi; | |
1191 | d_theta2 = iem->theta; | |
1192 | d_tpc2 = iem->tpc; | |
1193 | d_ntpc_ele2 = iem->ntpc_ele; | |
1194 | d_ntpc_pio2 = iem->ntpc_pio; | |
1195 | d_ntpc_kao2 = iem->ntpc_kao; | |
1196 | d_ntpc_pro2 = iem->ntpc_pro; | |
1197 | d_ntof_ele2 = iem->ntof_ele; | |
1198 | d_ntof_pio2 = iem->ntof_pio; | |
1199 | d_ntof_kao2 = iem->ntof_kao; | |
1200 | d_ntof_pro2 = iem->ntof_pro; | |
1201 | d_its2 = iem->its; | |
1202 | d_nits2 = iem->nits; | |
1203 | d_ntpc2 = iem->ntpc; | |
1204 | d_e2 = iem->e; | |
1205 | d_dphi2 = iem->dphi; | |
1206 | d_deta2 = iem->deta; | |
1207 | d_dcaxy2 = iem->dcaxy; | |
1208 | d_dcaz2 = iem->dcaz; | |
1209 | d_conv2 = iem->conv_flag; | |
1210 | ||
1211 | d_run = fkRunNumber; | |
1212 | // d_event = d_evt; | |
1213 | ||
1214 | ||
1215 | int icent = (int)(fkCentrality/10.0); | |
1216 | if(icent<0) icent=0; | |
1217 | if(icent>=NCENT) icent=NCENT-1; | |
1218 | if(pair_cut()==true){ | |
1219 | hmasspt[type][icent]->Fill(d_mass, d_ptpair); | |
1220 | hmasspt[type][10]->Fill(d_mass, d_ptpair); | |
1221 | } | |
1222 | ||
1223 | ////////////// pt weighting for mixed event | |
1224 | float weight = 1; | |
1225 | float k = 10000; | |
1226 | if(type==0 || type==1 || type==2){ | |
1227 | weight = 1; | |
1228 | }else if(type==3 || type==5 || type==6){ | |
1229 | k = 18.7*exp(-2.4*d_pt2)+1.24; | |
1230 | weight = 1+1/k; | |
1231 | }else if(type==4){ | |
1232 | k = 18.7*exp(-2.4*d_pt1)+1.24; | |
1233 | weight = 1+1/k; | |
1234 | } | |
1235 | if(pair_cut()==true){ | |
1236 | hmasspt_weight[type][icent]->Fill(d_mass, d_ptpair, weight); | |
1237 | hmasspt_weight[type][10]->Fill(d_mass, d_ptpair, weight); | |
1238 | } | |
1239 | ||
1240 | ||
1241 | //d_ntpair->Fill(); | |
1242 | ||
1243 | } | |
1244 | ||
1245 | //____________________________________________ | |
1246 | bool ana_sgl::PairTrackcut(vector<etrk>::iterator e1, vector<etrk>::iterator e2){ | |
1247 | /* | |
1248 | double p1 = sqrt(pow(e1.px,2)+ | |
1249 | pow(e1.py,2)+ | |
1250 | pow(e1.pz,2)); | |
1251 | ||
1252 | double p2 = sqrt(pow(e2.px,2)+ | |
1253 | pow(e2.py,2)+ | |
1254 | pow(e2.pz,2)); | |
1255 | ||
1256 | if( | |
1257 | (e1.e/p1>0.7 && e1.e/p1<1.3) || | |
1258 | (e2.e/p1>0.7 && e2.e/p2<1.3) | |
1259 | ){ | |
1260 | return true; | |
1261 | }else{ | |
1262 | return false; | |
1263 | } | |
1264 | o */ | |
1265 | return true; | |
1266 | ||
1267 | } | |
1268 | ||
1269 | //____________________________________________ | |
1270 | void ana_sgl::check_conversion_pairs(vector<etrk> &e1, vector<etrk> &e2){ | |
1271 | vector<etrk>::iterator iep; | |
1272 | vector<etrk>::iterator iem; | |
1273 | bool reject = false; | |
1274 | if(e1.size()>0 && e2.size()>0){ | |
1275 | for(iep = e1.begin(); iep != e1.end(); ++iep){ | |
1276 | reject = false; | |
1277 | for(iem = e2.begin(); iem != e2.end(); ++iem){ | |
1278 | double mass, phiv, px, py, pz, pt, e, phi, eta, cos, psi; | |
1279 | calc_vars(iep, iem, mass, phiv, px, py, pz, pt, e, phi, eta, cos, psi); | |
1280 | if(magnetic_field_mm==true){ //pp | |
1281 | if(phiv<0.6 && iep->phi-iem->phi<0){ // this depends on the magntic field | |
1282 | reject = true; | |
1283 | iem->conv_flag = 0; | |
1284 | } | |
1285 | }else{ | |
1286 | if(phiv>acos(-1.0)-0.6 && iep->phi-iem->phi>0){ // this depends on the magntic field | |
1287 | reject = true; | |
1288 | iem->conv_flag = 0; | |
1289 | } | |
1290 | } | |
1291 | } | |
1292 | if(reject==true) iep->conv_flag=0; | |
1293 | } | |
1294 | } | |
1295 | } | |
1296 | ||
1297 | //____________________________________________ | |
1298 | void ana_sgl::check_ghost_pairs(vector<etrk> &e1){ | |
1299 | vector<etrk>::iterator iep; | |
1300 | vector<etrk>::iterator iem; | |
1301 | bool reject = false; | |
1302 | if(e1.size()>1){ | |
1303 | for(iep = e1.begin(); iep != e1.end(); ++iep){ | |
1304 | reject = false; | |
1305 | vector<etrk>::iterator iep2=iep; | |
1306 | ++iep2; | |
1307 | for(iem = iep2; iem != e1.end(); ++iem){ | |
1308 | double mass, phiv, px, py, pz, pt, e, phi, eta, cos, psi; | |
1309 | calc_vars(iep, iem, mass, phiv, px, py, pz, pt, e, phi, eta, cos, psi); | |
1310 | if(mass<0.01){ | |
1311 | reject = true; | |
1312 | iem->ghost_flag = 0; | |
1313 | } | |
1314 | } | |
1315 | if(reject==true) iep->ghost_flag=0; | |
1316 | } | |
1317 | } | |
1318 | } | |
1319 | ||
1320 | //____________________________________________ | |
1321 | //void ana_sgl::calc_vars(etrk iep, etrk iem, double &mass, double &phiv, double &px, double &py, double&pz, | |
1322 | void ana_sgl::calc_vars(vector<etrk>::iterator iep, vector<etrk>::iterator iem, double &mass, double &phiv, double &px, double &py, double&pz, | |
1323 | double &pt, double &e, double &phi, double &eta, double &cos, double &psi){ | |
1324 | ||
1325 | px = iep->px+iem->px; | |
1326 | py = iep->py+iem->py; | |
1327 | pz = iep->pz+iem->pz; | |
1328 | pt = sqrt(px*px+py*py); | |
1329 | double d_ppair = sqrt(pt*pt+pz*pz); | |
1330 | static const double me=0.0005109989; | |
1331 | e = sqrt(me*me+iep->px*iep->px+iep->py*iep->py+iep->pz*iep->pz) | |
1332 | + sqrt(me*me+iem->px*iem->px+iem->py*iem->py+iem->pz*iem->pz); | |
1333 | ||
1334 | mass = e*e-px*px-py*py-pz*pz; | |
1335 | if(mass<0){ | |
1336 | mass = mass; | |
1337 | }else{ | |
1338 | mass = sqrt(mass); | |
1339 | } | |
1340 | ||
1341 | ||
1342 | ||
1343 | phi = atan2(py, px); | |
1344 | eta = -0.5*TMath::Log((d_ppair+pz)/(d_ppair-pz)); | |
1345 | double p1 = sqrt(pow(iep->px,2)+pow(iep->py,2)+pow(iep->pz,2)); | |
1346 | double p2 = sqrt(pow(iem->px,2)+pow(iem->py,2)+pow(iem->pz,2)); | |
1347 | cos = acos((iep->px*iem->px+iep->py*iem->py+iep->pz*iem->pz)/(p1*p2)); | |
1348 | ||
1349 | ||
1350 | double dtheta = iep->theta-iem->theta; | |
1351 | psi = asin(dtheta/cos); | |
1352 | ||
1353 | ||
1354 | //unit vector of (pep+pem) | |
1355 | float pl = d_ppair; | |
1356 | float ux = px/pl; | |
1357 | float uy = py/pl; | |
1358 | float uz = pz/pl; | |
1359 | float ax = uy/sqrt(ux*ux+uy*uy); | |
1360 | float ay = -ux/sqrt(ux*ux+uy*uy); | |
1361 | ||
1362 | //momentum of e+ and e- in (ax,ay,az) axis. Note that az=0 by | |
1363 | //definition. | |
1364 | float ptep = iep->px*ax + iep->py*ay; | |
1365 | float ptem = iem->px*ax + iem->py*ay; | |
1366 | ||
1367 | float pxep = iep->px; | |
1368 | float pyep = iep->py; | |
1369 | float pzep = iep->pz; | |
1370 | float pxem = iem->px; | |
1371 | float pyem = iem->py; | |
1372 | float pzem = iem->pz; | |
1373 | ||
1374 | ||
1375 | //vector product of pep X pem | |
1376 | float vpx = pyep*pzem - pzep*pyem; | |
1377 | float vpy = pzep*pxem - pxep*pzem; | |
1378 | float vpz = pxep*pyem - pyep*pxem; | |
1379 | float vp = sqrt(vpx*vpx+vpy*vpy+vpz*vpz); | |
1380 | float thev = acos(vpz/vp); | |
1381 | ||
1382 | //unit vector of pep X pem | |
1383 | float vx = vpx/vp; | |
1384 | float vy = vpy/vp; | |
1385 | float vz = vpz/vp; | |
1386 | ||
1387 | //The third axis defined by vector product (ux,uy,uz)X(vx,vy,vz) | |
1388 | float wx = uy*vz - uz*vy; | |
1389 | float wy = uz*vx - ux*vz; | |
1390 | float wz = ux*vy - uy*vx; | |
1391 | float wl = sqrt(wx*wx+wy*wy+wz*wz); | |
1392 | // by construction, (wx,wy,wz) must be a unit vector. | |
1393 | if(fabs(wl - 1.0) > 0.00001) cout << "Calculation error in W vector"<<endl; | |
1394 | // measure angle between (wx,wy,wz) and (ax,ay,0). The angle between them | |
1395 | // should be small if the pair is conversion | |
1396 | // | |
1397 | float cosPhiV = wx*ax + wy*ay; | |
1398 | phiv = acos(cosPhiV); | |
1399 | ||
1400 | } | |
1401 | ||
1402 | //____________________________________________ | |
1403 | bool ana_sgl::pair_cut(void){ | |
1404 | bool ret = true; | |
1405 | ||
1406 | ||
1407 | if(d_flag_phiv==true){ | |
1408 | if(d_run>169591 && d_phiv>acos(-1.0)-d_phiv_cut){ | |
1409 | ret = false; | |
1410 | } | |
1411 | if(d_run<169591 && d_phiv<d_phiv_cut){ | |
1412 | ret = false; | |
1413 | } | |
1414 | } | |
1415 | ||
1416 | ||
1417 | if(d_flag_emc_cut==true){ | |
1418 | if( !( | |
1419 | (d_e1/d_pt1>d_emc_low && d_e1/d_pt1<d_emc_high) || | |
1420 | (d_e2/d_pt2>d_emc_low && d_e2/d_pt2<d_emc_high) | |
1421 | ) | |
1422 | ){ | |
1423 | ret = false; | |
1424 | } | |
1425 | } | |
1426 | ||
1427 | ||
1428 | if(d_flag_pt_cut==true){ | |
1429 | if(!(d_pt1>d_pt_cut_low && d_pt2>d_pt_cut_low && | |
1430 | d_pt1<d_pt_cut_high && d_pt2<d_pt_cut_high) | |
1431 | ){ | |
1432 | ret = false; | |
1433 | } | |
1434 | } | |
1435 | ||
1436 | return ret; | |
1437 | ||
1438 | } | |
1439 | ||
1440 | //____________________________________________ | |
1441 | void ana_sgl::print_cuts(void){ | |
1442 | ||
1443 | cout<<" ********** lists of cuts ********** "<<endl; | |
1444 | cout<<d_tpc_dedx_low<<" < TPC dE/dx < "<<d_tpc_dedx_high<<endl; | |
1445 | ||
1446 | ||
1447 | if(d_flag_kaon_veto==true){ | |
1448 | cout<<d_dedx_kaon_veto_low<<" < TPC veto for Kaon < "<<d_dedx_kaon_veto_high<<endl; | |
1449 | }else{ | |
1450 | cout<<" No TPC Kaon Veto "<<endl; | |
1451 | } | |
1452 | ||
1453 | if(d_flag_proton_veto==true){ | |
1454 | cout<<d_dedx_proton_veto_low<<" < TPC veto for proton < "<<d_dedx_proton_veto_high<<endl; | |
1455 | }else{ | |
1456 | cout<<" No TPC Proton Veto "<<endl; | |
1457 | } | |
1458 | ||
1459 | if(d_flag_tof_cut==true){ | |
1460 | cout<<d_tof_low<<" < kTOFnSigmaEle < "<<d_tof_high<<endl; | |
1461 | }else{ | |
1462 | cout<<" No TOF cuts "<<endl; | |
1463 | } | |
1464 | if(d_flag_emc_cut==true){ | |
1465 | cout<<d_emc_low<<" < EMC E/p (or for pairs) < "<<d_emc_high<<endl; | |
1466 | }else{ | |
1467 | cout<<" No EMC cuts "<<endl; | |
1468 | } | |
1469 | if(d_flag_phiv==true){ | |
1470 | cout<<" Phiv cuts "<<d_phiv_cut<<endl; | |
1471 | }else{ | |
1472 | cout<<" No phiv cuts "<<endl; | |
1473 | } | |
1474 | ||
1475 | if(d_flag_pt_cut==true){ | |
1476 | cout<<d_pt_cut_low<<" < pT of electrons < "<<d_pt_cut_high<<endl; | |
1477 | }else{ | |
1478 | cout<<" No electron pt cuts "<<endl; | |
1479 | } | |
1480 | ||
1481 | if(d_conv_flag==true){ | |
1482 | cout<<" remove converson pairs from electron/positron lists "<<endl; | |
1483 | }else{ | |
1484 | cout<<" conversion candidates are in the pool."<<endl; | |
1485 | } | |
1486 | cout<<" ********** end of lists of cuts ********** "<<endl; | |
1487 | ||
1488 | if(sel_trigger==0){ | |
1489 | cout<<" MB "<<endl; | |
1490 | }else if(sel_trigger==1){ | |
1491 | cout<<" SemiCentral "<<endl; | |
1492 | }else if(sel_trigger==2){ | |
1493 | cout<<" Central "<<endl; | |
1494 | } | |
1495 | ||
1496 | } | |
1497 | ||
1498 |