]>
Commit | Line | Data |
---|---|---|
04a7657f | 1 | /* $Id$ */ |
2 | ||
3 | #include <AliPWG0Helper.h> | |
4 | ||
5 | #include <TParticle.h> | |
6 | #include <TParticlePDG.h> | |
847489f7 | 7 | #include <TH1.h> |
dd367a14 | 8 | #include <TH2.h> |
25db2d85 | 9 | #include <TH3.h> |
116083b1 | 10 | #include <TList.h> |
9cd015f7 | 11 | #include <TTree.h> |
12 | #include <TBranch.h> | |
13 | #include <TLeaf.h> | |
116083b1 | 14 | |
15 | #include <AliHeader.h> | |
16 | #include <AliStack.h> | |
04a7657f | 17 | |
18 | #include <AliLog.h> | |
19 | #include <AliESD.h> | |
770a1f1d | 20 | #include <AliESDEvent.h> |
04a7657f | 21 | #include <AliESDVertex.h> |
81be4ee8 | 22 | #include <AliESDRun.h> |
85737496 | 23 | #include <AliVertexerTracks.h> |
69b09e3b | 24 | #include <AliMultiplicity.h> |
04a7657f | 25 | |
116083b1 | 26 | #include <AliGenEventHeader.h> |
27 | #include <AliGenPythiaEventHeader.h> | |
28 | #include <AliGenCocktailEventHeader.h> | |
6b62a9c7 | 29 | #include <AliGenDPMjetEventHeader.h> |
69b09e3b | 30 | #include <AliESDVZERO.h> |
116083b1 | 31 | |
04a7657f | 32 | //____________________________________________________________________ |
33 | ClassImp(AliPWG0Helper) | |
34 | ||
a67484a6 | 35 | Int_t AliPWG0Helper::fgLastProcessType = -1; |
04a7657f | 36 | |
51f6de65 | 37 | //____________________________________________________________________ |
a6e0ebfe | 38 | Bool_t AliPWG0Helper::TestVertex(const AliESDVertex* vertex, AnalysisMode analysisMode, Bool_t debug) |
51f6de65 | 39 | { |
40 | // Checks if a vertex meets the needed quality criteria | |
41 | ||
42 | Float_t requiredZResolution = -1; | |
2701e5ae | 43 | if (analysisMode & kSPD || analysisMode & kTPCITS || analysisMode & kTPCSPD) |
51f6de65 | 44 | { |
d879a92d | 45 | // disable cut on resolution |
46 | requiredZResolution = 1000; | |
51f6de65 | 47 | } |
a7f69e56 | 48 | else if (analysisMode & kTPC) |
51f6de65 | 49 | requiredZResolution = 10.; |
50 | ||
51 | // check resolution | |
52 | Double_t zRes = vertex->GetZRes(); | |
53 | ||
54 | if (zRes > requiredZResolution) { | |
55 | if (debug) | |
56 | Printf("AliPWG0Helper::TestVertex: Resolution too poor %f (required: %f", zRes, requiredZResolution); | |
57 | return kFALSE; | |
58 | } | |
d879a92d | 59 | |
60 | if (vertex->IsFromVertexerZ()) | |
61 | { | |
62 | if (vertex->GetDispersion() > 0.02) | |
63 | { | |
64 | if (debug) | |
65 | Printf("AliPWG0Helper::TestVertex: Delta Phi too large in Vertexer Z: %f (required: %f", vertex->GetDispersion(), 0.02); | |
66 | return kFALSE; | |
67 | } | |
68 | } | |
51f6de65 | 69 | |
70 | return kTRUE; | |
71 | } | |
72 | ||
770a1f1d | 73 | //____________________________________________________________________ |
81be4ee8 | 74 | const AliESDVertex* AliPWG0Helper::GetVertex(const AliESDEvent* aEsd, AnalysisMode analysisMode, Bool_t debug) |
770a1f1d | 75 | { |
76 | // Get the vertex from the ESD and returns it if the vertex is valid | |
77 | // | |
78 | // Second argument decides which vertex is used (this selects | |
79 | // also the quality criteria that are applied) | |
80 | ||
81 | const AliESDVertex* vertex = 0; | |
81be4ee8 | 82 | if (analysisMode & kSPD) |
770a1f1d | 83 | { |
d1f50534 | 84 | vertex = aEsd->GetPrimaryVertexSPD(); |
ebf31fda | 85 | if (debug) |
86 | Printf("AliPWG0Helper::GetVertex: Returning SPD vertex"); | |
770a1f1d | 87 | } |
2701e5ae | 88 | else if (analysisMode & kTPCITS || analysisMode & kTPCSPD) |
81be4ee8 | 89 | { |
90 | vertex = aEsd->GetPrimaryVertexTracks(); | |
91 | if (debug) | |
92 | Printf("AliPWG0Helper::GetVertex: Returning vertex from tracks"); | |
68fa248f | 93 | if (!vertex || vertex->GetNContributors() <= 0) |
81be4ee8 | 94 | { |
95 | if (debug) | |
96 | Printf("AliPWG0Helper::GetVertex: Vertex from tracks has no contributors. Falling back to SPD vertex."); | |
97 | vertex = aEsd->GetPrimaryVertexSPD(); | |
98 | } | |
99 | } | |
a7f69e56 | 100 | else if (analysisMode & kTPC) |
770a1f1d | 101 | { |
d1f50534 | 102 | vertex = aEsd->GetPrimaryVertexTPC(); |
ebf31fda | 103 | if (debug) |
81be4ee8 | 104 | Printf("AliPWG0Helper::GetVertex: Returning vertex from TPC-only tracks"); |
770a1f1d | 105 | } |
106 | else | |
107 | Printf("AliPWG0Helper::GetVertex: ERROR: Invalid second argument %d", analysisMode); | |
108 | ||
ebf31fda | 109 | if (!vertex) { |
110 | if (debug) | |
111 | Printf("AliPWG0Helper::GetVertex: No vertex found in ESD"); | |
770a1f1d | 112 | return 0; |
ebf31fda | 113 | } |
770a1f1d | 114 | |
115 | // check Ncontributors | |
ebf31fda | 116 | if (vertex->GetNContributors() <= 0) { |
85737496 | 117 | if (debug){ |
118 | Printf("AliPWG0Helper::GetVertex: NContributors() <= 0: %d",vertex->GetNContributors()); | |
119 | Printf("AliPWG0Helper::GetVertex: NIndices(): %d",vertex->GetNIndices()); | |
51f6de65 | 120 | vertex->Print(); |
85737496 | 121 | } |
770a1f1d | 122 | return 0; |
ebf31fda | 123 | } |
770a1f1d | 124 | |
125 | // check resolution | |
126 | Double_t zRes = vertex->GetZRes(); | |
51f6de65 | 127 | if (zRes == 0) { |
128 | Printf("AliPWG0Helper::GetVertex: UNEXPECTED: resolution is 0."); | |
770a1f1d | 129 | return 0; |
ebf31fda | 130 | } |
131 | ||
132 | if (debug) | |
85737496 | 133 | { |
134 | Printf("AliPWG0Helper::GetVertex: Returning valid vertex: %s", vertex->GetTitle()); | |
135 | vertex->Print(); | |
136 | } | |
770a1f1d | 137 | |
138 | return vertex; | |
139 | } | |
140 | ||
04a7657f | 141 | //____________________________________________________________________ |
7584d357 | 142 | Bool_t AliPWG0Helper::IsPrimaryCharged(TParticle* aParticle, Int_t aTotalPrimaries, Bool_t adebug) |
04a7657f | 143 | { |
144 | // | |
25db2d85 | 145 | // this function checks if a particle from the event generator (i.e. among the nPrim particles in the stack) |
146 | // shall be counted as a primary particle | |
147 | // | |
04a7657f | 148 | // This function or a equivalent should be available in some common place of AliRoot |
149 | // | |
dd367a14 | 150 | // WARNING: Call this function only for particles that are among the particles from the event generator! |
151 | // --> stack->Particle(id) with id < stack->GetNprimary() | |
04a7657f | 152 | |
153 | // if the particle has a daughter primary, we do not want to count it | |
154 | if (aParticle->GetFirstDaughter() != -1 && aParticle->GetFirstDaughter() < aTotalPrimaries) | |
155 | { | |
7584d357 | 156 | if (adebug) |
25db2d85 | 157 | printf("Dropping particle because it has a daughter among the primaries.\n"); |
04a7657f | 158 | return kFALSE; |
159 | } | |
160 | ||
161 | Int_t pdgCode = TMath::Abs(aParticle->GetPdgCode()); | |
85737496 | 162 | |
04a7657f | 163 | |
164 | // skip quarks and gluon | |
165 | if (pdgCode <= 10 || pdgCode == 21) | |
166 | { | |
7584d357 | 167 | if (adebug) |
25db2d85 | 168 | printf("Dropping particle because it is a quark or gluon.\n"); |
04a7657f | 169 | return kFALSE; |
170 | } | |
171 | ||
85737496 | 172 | Int_t status = aParticle->GetStatusCode(); |
173 | // skip non final state particles.. | |
174 | if(status!=1){ | |
175 | if (adebug) | |
176 | printf("Dropping particle because it is not a final state particle.\n"); | |
177 | return kFALSE; | |
178 | } | |
179 | ||
04a7657f | 180 | if (strcmp(aParticle->GetName(),"XXX") == 0) |
181 | { | |
1c15d51a | 182 | Printf("WARNING: There is a particle named XXX (pdg code %d).", pdgCode); |
04a7657f | 183 | return kFALSE; |
184 | } | |
185 | ||
186 | TParticlePDG* pdgPart = aParticle->GetPDG(); | |
187 | ||
188 | if (strcmp(pdgPart->ParticleClass(),"Unknown") == 0) | |
189 | { | |
1c15d51a | 190 | Printf("WARNING: There is a particle with an unknown particle class (pdg code %d).", pdgCode); |
04a7657f | 191 | return kFALSE; |
192 | } | |
193 | ||
194 | if (pdgPart->Charge() == 0) | |
195 | { | |
7584d357 | 196 | if (adebug) |
25db2d85 | 197 | printf("Dropping particle because it is not charged.\n"); |
04a7657f | 198 | return kFALSE; |
199 | } | |
200 | ||
201 | return kTRUE; | |
202 | } | |
847489f7 | 203 | |
204 | //____________________________________________________________________ | |
29771dc8 | 205 | void AliPWG0Helper::CreateProjections(TH3* hist, Bool_t save) |
847489f7 | 206 | { |
207 | // create projections of 3d hists to all 2d combinations | |
208 | // the histograms are not returned, just use them from memory or use this to create them in a file | |
209 | ||
210 | TH1* proj = hist->Project3D("yx"); | |
211 | proj->SetXTitle(hist->GetXaxis()->GetTitle()); | |
212 | proj->SetYTitle(hist->GetYaxis()->GetTitle()); | |
29771dc8 | 213 | if (save) |
214 | proj->Write(); | |
847489f7 | 215 | |
216 | proj = hist->Project3D("zx"); | |
217 | proj->SetXTitle(hist->GetXaxis()->GetTitle()); | |
218 | proj->SetYTitle(hist->GetZaxis()->GetTitle()); | |
29771dc8 | 219 | if (save) |
220 | proj->Write(); | |
847489f7 | 221 | |
222 | proj = hist->Project3D("zy"); | |
223 | proj->SetXTitle(hist->GetYaxis()->GetTitle()); | |
224 | proj->SetYTitle(hist->GetZaxis()->GetTitle()); | |
29771dc8 | 225 | if (save) |
226 | proj->Write(); | |
847489f7 | 227 | } |
1afae8ff | 228 | |
229 | //____________________________________________________________________ | |
29771dc8 | 230 | void AliPWG0Helper::CreateDividedProjections(TH3* hist, TH3* hist2, const char* axis, Bool_t putErrors, Bool_t save) |
1afae8ff | 231 | { |
232 | // create projections of the 3d hists divides them | |
233 | // axis decides to which plane, if axis is 0 to all planes | |
234 | // the histograms are not returned, just use them from memory or use this to create them in a file | |
235 | ||
236 | if (axis == 0) | |
237 | { | |
29771dc8 | 238 | CreateDividedProjections(hist, hist2, "yx", putErrors, save); |
239 | CreateDividedProjections(hist, hist2, "zx", putErrors, save); | |
240 | CreateDividedProjections(hist, hist2, "zy", putErrors, save); | |
1afae8ff | 241 | |
242 | return; | |
243 | } | |
244 | ||
245 | TH1* proj = hist->Project3D(axis); | |
0ab29cfa | 246 | |
247 | if (strlen(axis) == 2) | |
248 | { | |
249 | proj->SetYTitle(GetAxisTitle(hist, axis[0])); | |
250 | proj->SetXTitle(GetAxisTitle(hist, axis[1])); | |
251 | } | |
252 | else if (strlen(axis) == 1) | |
253 | proj->SetXTitle(GetAxisTitle(hist, axis[0])); | |
1afae8ff | 254 | |
255 | TH1* proj2 = hist2->Project3D(axis); | |
0ab29cfa | 256 | if (strlen(axis) == 2) |
257 | { | |
258 | proj2->SetYTitle(GetAxisTitle(hist2, axis[0])); | |
259 | proj2->SetXTitle(GetAxisTitle(hist2, axis[1])); | |
260 | } | |
261 | else if (strlen(axis) == 1) | |
262 | proj2->SetXTitle(GetAxisTitle(hist2, axis[0])); | |
1afae8ff | 263 | |
264 | TH1* division = dynamic_cast<TH1*> (proj->Clone(Form("%s_div_%s", proj->GetName(), proj2->GetName()))); | |
29771dc8 | 265 | //printf("doing axis: %s, x axis has %d %d bins, min %f %f max %f %f\n", axis, division->GetNbinsX(), proj2->GetNbinsX(), division->GetXaxis()->GetBinLowEdge(1), proj2->GetXaxis()->GetBinLowEdge(1), division->GetXaxis()->GetBinUpEdge(division->GetNbinsX()), proj2->GetXaxis()->GetBinUpEdge(proj2->GetNbinsX())); |
266 | //printf("doing axis: %s, y axis has %d %d bins, min %f %f max %f %f\n", axis, division->GetNbinsY(), proj2->GetNbinsY(), division->GetYaxis()->GetBinLowEdge(1), proj2->GetYaxis()->GetBinLowEdge(1), division->GetYaxis()->GetBinUpEdge(division->GetNbinsY()), proj2->GetYaxis()->GetBinUpEdge(proj2->GetNbinsY())); | |
9e952c39 | 267 | division->Divide(proj, proj2, 1, 1, "B"); |
29771dc8 | 268 | division->SetTitle(Form("%s divided %s", proj->GetTitle(), proj2->GetTitle())); |
0ab29cfa | 269 | |
270 | if (putErrors) | |
271 | { | |
272 | division->Sumw2(); | |
273 | if (division->GetDimension() == 1) | |
274 | { | |
275 | Int_t nBins = division->GetNbinsX(); | |
29771dc8 | 276 | for (Int_t i = 1; i <= nBins; ++i) |
0ab29cfa | 277 | if (proj2->GetBinContent(i) != 0) |
278 | division->SetBinError(i, TMath::Sqrt(proj->GetBinContent(i)) / proj2->GetBinContent(i)); | |
279 | } | |
280 | else if (division->GetDimension() == 2) | |
281 | { | |
282 | Int_t nBinsX = division->GetNbinsX(); | |
283 | Int_t nBinsY = division->GetNbinsY(); | |
29771dc8 | 284 | for (Int_t i = 1; i <= nBinsX; ++i) |
285 | for (Int_t j = 1; j <= nBinsY; ++j) | |
0ab29cfa | 286 | if (proj2->GetBinContent(i, j) != 0) |
287 | division->SetBinError(i, j, TMath::Sqrt(proj->GetBinContent(i, j)) / proj2->GetBinContent(i, j)); | |
288 | } | |
289 | } | |
29771dc8 | 290 | |
291 | if (save) | |
292 | { | |
293 | proj->Write(); | |
294 | proj2->Write(); | |
295 | division->Write(); | |
296 | } | |
1afae8ff | 297 | } |
298 | ||
92d2d8ad | 299 | //____________________________________________________________________ |
25db2d85 | 300 | const char* AliPWG0Helper::GetAxisTitle(TH3* hist, const char axis) |
92d2d8ad | 301 | { |
302 | // returns the title of the axis given in axis (x, y, z) | |
303 | ||
304 | if (axis == 'x') | |
305 | return hist->GetXaxis()->GetTitle(); | |
306 | else if (axis == 'y') | |
307 | return hist->GetYaxis()->GetTitle(); | |
308 | else if (axis == 'z') | |
309 | return hist->GetZaxis()->GetTitle(); | |
310 | ||
311 | return 0; | |
312 | } | |
116083b1 | 313 | |
6b62a9c7 | 314 | |
1c15d51a | 315 | AliPWG0Helper::MCProcessType AliPWG0Helper::GetPythiaEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) { |
6b62a9c7 | 316 | |
317 | AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader); | |
318 | ||
319 | if (!pythiaGenHeader) { | |
320 | printf("AliPWG0Helper::GetProcessType : Unknown gen Header type). \n"); | |
321 | return kInvalidProcess; | |
322 | } | |
323 | ||
324 | ||
325 | Int_t pythiaType = pythiaGenHeader->ProcessType(); | |
a67484a6 | 326 | fgLastProcessType = pythiaType; |
6b62a9c7 | 327 | MCProcessType globalType = kInvalidProcess; |
328 | ||
329 | ||
330 | if (adebug) { | |
331 | printf("AliPWG0Helper::GetProcessType : Pythia process type found: %d \n",pythiaType); | |
332 | } | |
333 | ||
334 | ||
335 | if(pythiaType==92||pythiaType==93){ | |
336 | globalType = kSD; | |
337 | } | |
338 | else if(pythiaType==94){ | |
339 | globalType = kDD; | |
340 | } | |
341 | //else if(pythiaType != 91){ // also exclude elastic to be sure... CKB?? | |
342 | else { | |
343 | globalType = kND; | |
344 | } | |
345 | return globalType; | |
346 | } | |
347 | ||
348 | ||
1c15d51a | 349 | AliPWG0Helper::MCProcessType AliPWG0Helper::GetDPMjetEventProcessType(AliGenEventHeader* aHeader, Bool_t adebug) { |
116083b1 | 350 | // |
351 | // get the process type of the event. | |
352 | // | |
353 | ||
354 | // can only read pythia headers, either directly or from cocktalil header | |
6b62a9c7 | 355 | AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(aHeader); |
356 | ||
357 | if (!dpmJetGenHeader) { | |
358 | printf("AliPWG0Helper::GetDPMjetProcessType : Unknown header type (not DPMjet or). \n"); | |
359 | return kInvalidProcess; | |
360 | } | |
361 | ||
362 | Int_t dpmJetType = dpmJetGenHeader->ProcessType(); | |
a67484a6 | 363 | fgLastProcessType = dpmJetType; |
6b62a9c7 | 364 | MCProcessType globalType = kInvalidProcess; |
365 | ||
366 | ||
367 | if (adebug) { | |
368 | printf("AliPWG0Helper::GetDPMJetProcessType : DPMJet process type found: %d \n",dpmJetType); | |
369 | } | |
370 | ||
371 | ||
81be4ee8 | 372 | if (dpmJetType == 1 || dpmJetType == 4) { // explicitly inelastic plus central diffraction |
6b62a9c7 | 373 | globalType = kND; |
374 | } | |
81be4ee8 | 375 | else if (dpmJetType==5 || dpmJetType==6) { |
6b62a9c7 | 376 | globalType = kSD; |
377 | } | |
81be4ee8 | 378 | else if (dpmJetType==7) { |
6b62a9c7 | 379 | globalType = kDD; |
380 | } | |
381 | return globalType; | |
382 | } | |
383 | ||
384 | ||
1c15d51a | 385 | AliPWG0Helper::MCProcessType AliPWG0Helper::GetEventProcessType(AliHeader* aHeader, Bool_t adebug) { |
6b62a9c7 | 386 | // |
387 | // get the process type of the event. | |
388 | // | |
389 | ||
390 | ||
391 | // Check for simple headers first | |
392 | ||
116083b1 | 393 | AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(aHeader->GenEventHeader()); |
6b62a9c7 | 394 | if (pythiaGenHeader) { |
395 | return GetPythiaEventProcessType(pythiaGenHeader,adebug); | |
396 | } | |
116083b1 | 397 | |
6b62a9c7 | 398 | AliGenDPMjetEventHeader* dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(aHeader->GenEventHeader()); |
399 | if (dpmJetGenHeader) { | |
400 | return GetDPMjetEventProcessType(dpmJetGenHeader,adebug); | |
401 | } | |
402 | ||
116083b1 | 403 | |
6b62a9c7 | 404 | // check for cocktail |
116083b1 | 405 | |
6b62a9c7 | 406 | AliGenCocktailEventHeader* genCocktailHeader = dynamic_cast<AliGenCocktailEventHeader*>(aHeader->GenEventHeader()); |
407 | if (!genCocktailHeader) { | |
408 | printf("AliPWG0Helper::GetProcessType : Unknown header type (not Pythia or Cocktail). \n"); | |
409 | return kInvalidProcess; | |
410 | } | |
116083b1 | 411 | |
6b62a9c7 | 412 | TList* headerList = genCocktailHeader->GetHeaders(); |
413 | if (!headerList) { | |
414 | return kInvalidProcess; | |
415 | } | |
416 | ||
417 | for (Int_t i=0; i<headerList->GetEntries(); i++) { | |
418 | ||
419 | pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(headerList->At(i)); | |
420 | if (pythiaGenHeader) { | |
421 | return GetPythiaEventProcessType(pythiaGenHeader,adebug); | |
116083b1 | 422 | } |
423 | ||
6b62a9c7 | 424 | dpmJetGenHeader = dynamic_cast<AliGenDPMjetEventHeader*>(headerList->At(i)); |
425 | if (dpmJetGenHeader) { | |
426 | return GetDPMjetEventProcessType(dpmJetGenHeader,adebug); | |
116083b1 | 427 | } |
428 | } | |
6b62a9c7 | 429 | return kInvalidProcess; |
430 | } | |
116083b1 | 431 | |
116083b1 | 432 | |
116083b1 | 433 | |
434 | //____________________________________________________________________ | |
435 | TParticle* AliPWG0Helper::FindPrimaryMother(AliStack* stack, Int_t label) | |
436 | { | |
437 | // | |
438 | // Finds the first mother among the primary particles of the particle identified by <label>, | |
439 | // i.e. the primary that "caused" this particle | |
440 | // | |
441 | ||
442 | Int_t motherLabel = FindPrimaryMotherLabel(stack, label); | |
443 | if (motherLabel < 0) | |
444 | return 0; | |
445 | ||
446 | return stack->Particle(motherLabel); | |
447 | } | |
448 | ||
449 | //____________________________________________________________________ | |
450 | Int_t AliPWG0Helper::FindPrimaryMotherLabel(AliStack* stack, Int_t label) | |
451 | { | |
452 | // | |
453 | // Finds the first mother among the primary particles of the particle identified by <label>, | |
454 | // i.e. the primary that "caused" this particle | |
455 | // | |
456 | // returns its label | |
457 | // | |
458 | ||
459 | Int_t nPrim = stack->GetNprimary(); | |
460 | ||
461 | while (label >= nPrim) | |
462 | { | |
463 | //printf("Particle %d (pdg %d) is not a primary. Let's check its mother %d\n", label, mother->GetPdgCode(), mother->GetMother(0)); | |
464 | ||
465 | TParticle* particle = stack->Particle(label); | |
466 | if (!particle) | |
467 | { | |
468 | AliDebugGeneral("FindPrimaryMother", AliLog::kError, Form("UNEXPECTED: particle with label %d not found in stack.", label)); | |
469 | return -1; | |
470 | } | |
471 | ||
472 | // find mother | |
473 | if (particle->GetMother(0) < 0) | |
474 | { | |
475 | AliDebugGeneral("FindPrimaryMother", AliLog::kError, Form("UNEXPECTED: Could not find mother of secondary particle %d.", label)); | |
476 | return -1; | |
477 | } | |
478 | ||
479 | label = particle->GetMother(0); | |
480 | } | |
481 | ||
482 | return label; | |
483 | } | |
9cd015f7 | 484 | |
dd367a14 | 485 | //____________________________________________________________________ |
486 | void AliPWG0Helper::NormalizeToBinWidth(TH1* hist) | |
487 | { | |
488 | // | |
489 | // normalizes a 1-d histogram to its bin width | |
490 | // | |
491 | ||
492 | for (Int_t i=1; i<=hist->GetNbinsX(); ++i) | |
493 | { | |
494 | hist->SetBinContent(i, hist->GetBinContent(i) / hist->GetBinWidth(i)); | |
495 | hist->SetBinError(i, hist->GetBinError(i) / hist->GetBinWidth(i)); | |
496 | } | |
497 | } | |
498 | ||
499 | //____________________________________________________________________ | |
500 | void AliPWG0Helper::NormalizeToBinWidth(TH2* hist) | |
501 | { | |
502 | // | |
503 | // normalizes a 2-d histogram to its bin width (x width * y width) | |
504 | // | |
505 | ||
506 | for (Int_t i=1; i<=hist->GetNbinsX(); ++i) | |
507 | for (Int_t j=1; j<=hist->GetNbinsY(); ++j) | |
508 | { | |
509 | Double_t factor = hist->GetXaxis()->GetBinWidth(i) * hist->GetYaxis()->GetBinWidth(j); | |
510 | hist->SetBinContent(i, j, hist->GetBinContent(i, j) / factor); | |
511 | hist->SetBinError(i, j, hist->GetBinError(i, j) / factor); | |
512 | } | |
513 | } | |
d1f50534 | 514 | |
ff8c4f30 | 515 | //____________________________________________________________________ |
81be4ee8 | 516 | void AliPWG0Helper::PrintConf(AnalysisMode analysisMode, AliTriggerAnalysis::Trigger trigger, AliPWG0Helper::DiffTreatment diffTreatment) |
d1f50534 | 517 | { |
518 | // | |
519 | // Prints the given configuration | |
520 | // | |
521 | ||
a7f69e56 | 522 | TString str(">>>> Running with >"); |
523 | ||
524 | if (analysisMode & kSPD) | |
525 | str += "SPD-only"; | |
526 | ||
d879a92d | 527 | if (analysisMode & kSPDOnlyL0) |
528 | str += " (only L0 clusters)"; | |
529 | ||
a7f69e56 | 530 | if (analysisMode & kTPC) |
531 | str += "TPC-only"; | |
532 | ||
533 | if (analysisMode & kTPCITS) | |
534 | str += "Global tracking"; | |
2701e5ae | 535 | |
536 | if (analysisMode & kTPCSPD) | |
537 | str += "Tracks and tracklets"; | |
d1f50534 | 538 | |
a7f69e56 | 539 | if (analysisMode & kFieldOn) |
d1f50534 | 540 | { |
a7f69e56 | 541 | str += " (with field)"; |
d1f50534 | 542 | } |
a7f69e56 | 543 | else |
544 | str += " (WITHOUT field)"; | |
545 | ||
546 | str += "< and trigger >"; | |
70fdd197 | 547 | str += AliTriggerAnalysis::GetTriggerName(trigger); |
81be4ee8 | 548 | str += "< and diffractive treatment >"; |
549 | ||
550 | switch (diffTreatment) | |
551 | { | |
552 | case kMCFlags: | |
553 | str += "MC flags"; | |
554 | break; | |
555 | ||
556 | case kUA5Cuts: | |
557 | str += "UA5 cuts"; | |
558 | break; | |
559 | ||
560 | case kE710Cuts: | |
561 | str += "E710 cuts"; | |
562 | break; | |
563 | ||
564 | case kALICEHadronLevel: | |
565 | str += "ALICE hadron level"; | |
566 | break; | |
567 | } | |
568 | ||
a7f69e56 | 569 | str += "< <<<<"; |
d1f50534 | 570 | |
571 | Printf("%s", str.Data()); | |
572 | } | |
573 | ||
81be4ee8 | 574 | //____________________________________________________________________ |
575 | Double_t AliPWG0Helper::Rapidity(Double_t pt, Double_t pz, Double_t m) | |
576 | { | |
577 | // | |
578 | // calculates rapidity keeping the sign in case E == pz | |
579 | // | |
580 | ||
581 | Double_t energy = TMath::Sqrt(pt*pt+pz*pz+m*m); | |
582 | if (energy != TMath::Abs(pz)) | |
583 | return 0.5*TMath::Log((energy+pz)/(energy-pz)); | |
584 | ||
585 | Printf("W- mt=0"); | |
586 | return TMath::Sign(1.e30,pz); | |
587 | } | |
588 | ||
589 | //____________________________________________________________________ | |
590 | Bool_t AliPWG0Helper::IsHadronLevelSingleDiffractive(AliStack* stack, Float_t cms, Float_t xiMin, Float_t xiMax) | |
591 | { | |
592 | // | |
593 | // return if process is single diffractive on hadron level | |
594 | // | |
595 | // xiMax and xiMin cut on M^2/s | |
596 | // | |
597 | // Based on code from Martin Poghoysan | |
598 | // | |
599 | ||
600 | TParticle* part1 = 0; | |
601 | TParticle* part2 = 0; | |
602 | ||
603 | Double_t smallestY = 1e10; | |
604 | Double_t largestY = -1e10; | |
605 | ||
606 | for (Int_t iParticle = 0; iParticle < stack->GetNprimary(); iParticle++) | |
607 | { | |
608 | TParticle* part = stack->Particle(iParticle); | |
609 | if (!part) | |
610 | continue; | |
611 | ||
612 | Int_t pdg = TMath::Abs(part->GetPdgCode()); | |
613 | ||
614 | Int_t child1 = part->GetFirstDaughter(); | |
615 | Int_t ist = part->GetStatusCode(); | |
616 | ||
617 | Int_t mfl = Int_t (pdg / TMath::Power(10, Int_t(TMath::Log10(pdg)))); | |
618 | if (child1 > -1 || ist != 1) | |
619 | mfl = 0; // select final state charm and beauty | |
620 | if (!(stack->IsPhysicalPrimary(iParticle) || pdg == 111 || pdg == 3212 || pdg==3124 || mfl >= 4)) | |
621 | continue; | |
622 | Int_t imother = part->GetFirstMother(); | |
623 | if (imother>0) | |
624 | { | |
625 | TParticle *partM = stack->Particle(imother); | |
626 | Int_t pdgM=TMath::Abs(partM->GetPdgCode()); | |
627 | if (pdgM==111 || pdgM==3124 || pdgM==3212) | |
628 | continue; | |
629 | } | |
630 | ||
631 | Double_t y = 0; | |
632 | ||
633 | // fix for problem with getting mass of particle 3124 | |
634 | if (pdg != 3124) | |
635 | y = Rapidity(part->Pt(), part->Pz(), part->GetMass()); | |
636 | else | |
637 | y = Rapidity(part->Pt(), part->Pz(), 1.5195); | |
638 | ||
639 | if (y < smallestY) | |
640 | { | |
641 | smallestY = y; | |
642 | part1 = part; | |
643 | } | |
644 | ||
645 | if (y > largestY) | |
646 | { | |
647 | largestY = y; | |
648 | part2 = part; | |
649 | } | |
650 | } | |
651 | ||
652 | if (part1 == 0 || part2 == 0) | |
653 | return kFALSE; | |
654 | ||
655 | Int_t pdg1 = part1->GetPdgCode(); | |
656 | Int_t pdg2 = part2->GetPdgCode(); | |
657 | ||
658 | Double_t pt1 = part1->Pt(); | |
659 | Double_t pt2 = part2->Pt(); | |
660 | Double_t pz1 = part1->Pz(); | |
661 | Double_t pz2 = part2->Pz(); | |
662 | ||
663 | Double_t y1 = TMath::Abs(Rapidity(pt1, pz1, 0.938)); | |
664 | Double_t y2 = TMath::Abs(Rapidity(pt2, pz2, 0.938)); | |
665 | ||
666 | Int_t arm = -99999; | |
667 | if (pdg1 == 2212 && pdg2 == 2212) | |
668 | { | |
669 | if (y1 > y2) | |
670 | arm = 0; | |
671 | else | |
672 | arm = 1; | |
673 | } | |
674 | else if (pdg1 == 2212) | |
675 | arm = 0; | |
676 | else if (pdg2 == 2212) | |
677 | arm = 1; | |
678 | ||
679 | Double_t M02s = 1. - 2 * part1->Energy() / cms; | |
680 | Double_t M12s = 1. - 2 * part2->Energy() / cms; | |
681 | ||
682 | if (arm == 0 && M02s > xiMin && M02s < xiMax) | |
683 | return kTRUE; | |
684 | else if (arm == 1 && M12s > xiMin && M12s < xiMax) | |
685 | return kTRUE; | |
686 | ||
687 | return kFALSE; | |
688 | } | |
689 | ||
690 | //____________________________________________________________________ | |
691 | AliPWG0Helper::MCProcessType AliPWG0Helper::GetEventProcessType(AliESDEvent* esd, AliHeader* header, AliStack* stack, AliPWG0Helper::DiffTreatment diffTreatment) | |
692 | { | |
693 | // | |
694 | // get process type | |
695 | // | |
696 | // diffTreatment: | |
697 | // kMCFlags: use MC flags | |
698 | // kUA5Cuts: SD events are those that have the MC SD flag and fulfill M^2/s < 0.05; DD from MC flag; Remainder is ND | |
699 | // kE710Cuts: SD events are those that have the MC SD flag and fulfill 2 < M^2 < 0.05s; DD from MC flag; Remainder is ND | |
700 | // kALICEHadronLevel: SD events are those that fulfill M^2/s < 0.05; DD from MC flag; Remainder is ND | |
701 | // | |
702 | ||
703 | MCProcessType mcProcessType = GetEventProcessType(header); | |
704 | ||
705 | if (diffTreatment == kMCFlags) | |
706 | return mcProcessType; | |
707 | ||
708 | Float_t cms = esd->GetESDRun()->GetBeamEnergy(); | |
709 | if (esd->GetESDRun()->IsBeamEnergyIsSqrtSHalfGeV()) | |
710 | cms *= 2; | |
711 | //Printf("cms = %f", cms); | |
712 | ||
713 | if (diffTreatment == kUA5Cuts && mcProcessType == kSD) | |
714 | { | |
715 | if (IsHadronLevelSingleDiffractive(stack, cms, 0, 0.05)) | |
716 | return kSD; | |
717 | } | |
718 | else if (diffTreatment == kE710Cuts && mcProcessType == kSD) | |
719 | { | |
720 | if (IsHadronLevelSingleDiffractive(stack, cms, 2. / cms / cms, 0.05)) | |
721 | return kSD; | |
722 | } | |
723 | else if (diffTreatment == kALICEHadronLevel) | |
724 | { | |
725 | if (IsHadronLevelSingleDiffractive(stack, cms, 0, 0.05)) | |
726 | return kSD; | |
727 | } | |
728 | ||
729 | if (mcProcessType == kSD) | |
730 | return kND; | |
731 | ||
732 | return mcProcessType; | |
733 | } |