]>
Commit | Line | Data |
---|---|---|
7cd4e982 | 1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * | |
3 | * * | |
4 | * Author: The ALICE Off-line Project. * | |
5 | * Contributors are mentioned in the code where appropriate. * | |
6 | * * | |
7 | * Permission to use, copy, modify and distribute this software and its * | |
8 | * documentation strictly for non-commercial purposes is hereby granted * | |
9 | * without fee, provided that the above copyright notice appears in all * | |
10 | * copies and that both the copyright notice and this permission notice * | |
11 | * appear in the supporting documentation. The authors make no claims * | |
12 | * about the suitability of this software for any purpose. It is * | |
13 | * provided "as is" without express or implied warranty. * | |
14 | **************************************************************************/ | |
7cd4e982 | 15 | |
16 | //_________________________________________________________________________ | |
17 | // Class for analysis utils for MC data | |
18 | // stored in stack or event header. | |
19 | // Contains: | |
20 | // - method to check the origin of a given track/cluster | |
21 | // - method to obtain the generated jets | |
22 | // | |
23 | //*-- Author: Gustavo Conesa (LNF-INFN) | |
24 | ////////////////////////////////////////////////////////////////////////////// | |
f8006433 | 25 | |
7cd4e982 | 26 | |
27 | // --- ROOT system --- | |
28 | #include <TMath.h> | |
29 | #include <TList.h> | |
30 | #include "TParticle.h" | |
8dacfd76 | 31 | #include "TDatabasePDG.h" |
c4a7d28a | 32 | #include "TVector3.h" |
7cd4e982 | 33 | |
34 | //---- ANALYSIS system ---- | |
35 | #include "AliMCAnalysisUtils.h" | |
36 | #include "AliCaloTrackReader.h" | |
37 | #include "AliStack.h" | |
38 | #include "AliGenPythiaEventHeader.h" | |
39 | #include "AliAODMCParticle.h" | |
40 | ||
f8006433 | 41 | ClassImp(AliMCAnalysisUtils) |
7cd4e982 | 42 | |
c5693f62 | 43 | //________________________________________ |
44 | AliMCAnalysisUtils::AliMCAnalysisUtils() : | |
45 | TObject(), | |
46 | fCurrentEvent(-1), | |
47 | fDebug(-1), | |
48 | fJetsList(new TList), | |
0d360687 | 49 | fMCGenerator(kPythia), |
50 | fMCGeneratorString("PYTHIA"), | |
51 | fDaughMom(), fDaughMom2(), | |
52 | fMotherMom(), fGMotherMom() | |
7cd4e982 | 53 | { |
54 | //Ctor | |
55 | } | |
7cd4e982 | 56 | |
c5693f62 | 57 | //_______________________________________ |
7cd4e982 | 58 | AliMCAnalysisUtils::~AliMCAnalysisUtils() |
59 | { | |
60 | // Remove all pointers. | |
61 | ||
0d360687 | 62 | if (fJetsList) |
63 | { | |
7cd4e982 | 64 | fJetsList->Clear(); |
65 | delete fJetsList ; | |
66 | } | |
67 | } | |
68 | ||
c5693f62 | 69 | //_____________________________________________________________________________________________ |
8a2dbbff | 70 | Int_t AliMCAnalysisUtils::CheckCommonAncestor(Int_t index1, Int_t index2, |
c5693f62 | 71 | const AliCaloTrackReader* reader, |
72 | Int_t & ancPDG, Int_t & ancStatus, | |
73 | TLorentzVector & momentum, TVector3 & prodVertex) | |
840124fe | 74 | { |
0d360687 | 75 | // Check the first common ancestor of 2 clusters, given the most likely labels of the primaries generating such clusters. |
76 | ||
48a0baf4 | 77 | Int_t label1[100]; |
78 | Int_t label2[100]; | |
79 | label1[0]= index1; | |
80 | label2[0]= index2; | |
81 | Int_t counter1 = 0; | |
82 | Int_t counter2 = 0; | |
840124fe | 83 | |
2644ead9 | 84 | if(label1[0]==label2[0]) |
85 | { | |
48a0baf4 | 86 | //printf("AliMCAnalysisUtils::CheckCommonAncestor() - Already the same label: %d\n",label1[0]); |
87 | counter1=1; | |
88 | counter2=1; | |
89 | } | |
2644ead9 | 90 | else |
91 | { | |
92 | if(reader->ReadAODMCParticles()) | |
93 | { | |
94 | TClonesArray * mcparticles = reader->GetAODMCParticles(); | |
48a0baf4 | 95 | |
96 | Int_t label=label1[0]; | |
2644ead9 | 97 | if(label < 0) return -1; |
98 | ||
36769d30 | 99 | while(label > -1 && counter1 < 99) |
100 | { | |
48a0baf4 | 101 | counter1++; |
102 | AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label); | |
36769d30 | 103 | if(mom) |
104 | { | |
840124fe | 105 | label = mom->GetMother() ; |
106 | label1[counter1]=label; | |
48a0baf4 | 107 | } |
108 | //printf("\t counter %d, label %d\n", counter1,label); | |
109 | } | |
36769d30 | 110 | |
48a0baf4 | 111 | //printf("Org label2=%d,\n",label2[0]); |
112 | label=label2[0]; | |
2644ead9 | 113 | if(label < 0) return -1; |
36769d30 | 114 | |
115 | while(label > -1 && counter2 < 99) | |
116 | { | |
48a0baf4 | 117 | counter2++; |
118 | AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label); | |
36769d30 | 119 | if(mom) |
120 | { | |
48a0baf4 | 121 | label = mom->GetMother() ; |
122 | label2[counter2]=label; | |
123 | } | |
124 | //printf("\t counter %d, label %d\n", counter2,label); | |
125 | } | |
126 | }//AOD MC | |
36769d30 | 127 | else |
128 | { //Kine stack from ESDs | |
48a0baf4 | 129 | AliStack * stack = reader->GetStack(); |
36769d30 | 130 | |
48a0baf4 | 131 | Int_t label=label1[0]; |
36769d30 | 132 | while(label > -1 && counter1 < 99) |
133 | { | |
48a0baf4 | 134 | counter1++; |
135 | TParticle * mom = stack->Particle(label); | |
36769d30 | 136 | if(mom) |
137 | { | |
48a0baf4 | 138 | label = mom->GetFirstMother() ; |
139 | label1[counter1]=label; | |
140 | } | |
141 | //printf("\t counter %d, label %d\n", counter1,label); | |
142 | } | |
36769d30 | 143 | |
48a0baf4 | 144 | //printf("Org label2=%d,\n",label2[0]); |
36769d30 | 145 | |
48a0baf4 | 146 | label=label2[0]; |
36769d30 | 147 | while(label > -1 && counter2 < 99) |
148 | { | |
48a0baf4 | 149 | counter2++; |
150 | TParticle * mom = stack->Particle(label); | |
36769d30 | 151 | if(mom) |
152 | { | |
48a0baf4 | 153 | label = mom->GetFirstMother() ; |
154 | label2[counter2]=label; | |
155 | } | |
156 | //printf("\t counter %d, label %d\n", counter2,label); | |
157 | } | |
158 | }// Kine stack from ESDs | |
159 | }//First labels not the same | |
160 | ||
36769d30 | 161 | if((counter1==99 || counter2==99) && fDebug >=0) |
162 | printf("AliMCAnalysisUtils::CheckCommonAncestor() - Genealogy too large c1: %d, c2= %d\n", counter1, counter2); | |
48a0baf4 | 163 | //printf("CheckAncestor:\n"); |
36769d30 | 164 | |
48a0baf4 | 165 | Int_t commonparents = 0; |
166 | Int_t ancLabel = -1; | |
167 | //printf("counters %d %d \n",counter1, counter2); | |
36769d30 | 168 | for (Int_t c1 = 0; c1 < counter1; c1++) |
169 | { | |
170 | for (Int_t c2 = 0; c2 < counter2; c2++) | |
171 | { | |
172 | if(label1[c1]==label2[c2] && label1[c1]>-1) | |
173 | { | |
48a0baf4 | 174 | ancLabel = label1[c1]; |
175 | commonparents++; | |
36769d30 | 176 | |
177 | if(reader->ReadAODMCParticles()) | |
178 | { | |
2644ead9 | 179 | AliAODMCParticle * mom = (AliAODMCParticle *) reader->GetAODMCParticles()->At(label1[c1]); |
36769d30 | 180 | |
181 | if (mom) | |
182 | { | |
48a0baf4 | 183 | ancPDG = mom->GetPdgCode(); |
184 | ancStatus = mom->GetStatus(); | |
952615e5 | 185 | momentum.SetPxPyPzE(mom->Px(),mom->Py(),mom->Pz(),mom->E()); |
c4a7d28a | 186 | prodVertex.SetXYZ(mom->Xv(),mom->Yv(),mom->Zv()); |
48a0baf4 | 187 | } |
188 | } | |
36769d30 | 189 | else |
190 | { | |
48a0baf4 | 191 | TParticle * mom = (reader->GetStack())->Particle(label1[c1]); |
36769d30 | 192 | |
193 | if (mom) | |
194 | { | |
48a0baf4 | 195 | ancPDG = mom->GetPdgCode(); |
196 | ancStatus = mom->GetStatusCode(); | |
952615e5 | 197 | mom->Momentum(momentum); |
c4a7d28a | 198 | prodVertex.SetXYZ(mom->Vx(),mom->Vy(),mom->Vz()); |
48a0baf4 | 199 | } |
200 | } | |
36769d30 | 201 | |
48a0baf4 | 202 | //First ancestor found, end the loops |
203 | counter1=0; | |
204 | counter2=0; | |
205 | }//Ancestor found | |
206 | }//second cluster loop | |
207 | }//first cluster loop | |
208 | ||
36769d30 | 209 | if(ancLabel < 0) |
210 | { | |
211 | ancPDG = -10000; | |
212 | ancStatus = -10000; | |
213 | momentum.SetXYZT(0,0,0,0); | |
214 | prodVertex.SetXYZ(-10,-10,-10); | |
215 | } | |
216 | ||
48a0baf4 | 217 | return ancLabel; |
218 | } | |
219 | ||
3e57baab | 220 | //________________________________________________________________________________________ |
8a2dbbff | 221 | Int_t AliMCAnalysisUtils::CheckOrigin(const Int_t * label, Int_t nlabels, |
3e57baab | 222 | const AliCaloTrackReader* reader, TString calorimeter) |
840124fe | 223 | { |
0d360687 | 224 | // Play with the montecarlo particles if available. |
225 | ||
3c75ddf7 | 226 | Int_t tag = 0; |
227 | ||
840124fe | 228 | if(nlabels<=0) { |
229 | printf("AliMCAnalysisUtils::CheckOrigin(nlabel<=0) - No MC labels available, please check!!!\n"); | |
230 | return kMCBadLabel; | |
231 | } | |
232 | ||
3e57baab | 233 | TObjArray* arrayCluster = 0; |
234 | if ( calorimeter == "EMCAL" ) arrayCluster = reader->GetEMCALClusters(); | |
235 | else if ( calorimeter == "PHOS" ) arrayCluster= reader->GetPHOSClusters(); | |
236 | ||
3c75ddf7 | 237 | //Select where the information is, ESD-galice stack or AOD mcparticles branch |
238 | if(reader->ReadStack()){ | |
3e57baab | 239 | tag = CheckOriginInStack(label, nlabels, reader->GetStack(), arrayCluster); |
3c75ddf7 | 240 | } |
241 | else if(reader->ReadAODMCParticles()){ | |
3e57baab | 242 | tag = CheckOriginInAOD(label, nlabels, reader->GetAODMCParticles(),arrayCluster); |
3c75ddf7 | 243 | } |
244 | ||
245 | return tag ; | |
840124fe | 246 | } |
247 | ||
3e57baab | 248 | //____________________________________________________________________________________________________ |
249 | Int_t AliMCAnalysisUtils::CheckOrigin(Int_t label, const AliCaloTrackReader* reader, TString calorimeter) | |
840124fe | 250 | { |
0d360687 | 251 | // Play with the montecarlo particles if available. |
252 | ||
3c75ddf7 | 253 | Int_t tag = 0; |
b5ef1905 | 254 | |
255 | if(label<0) { | |
898c9d44 | 256 | printf("AliMCAnalysisUtils::CheckOrigin(label<0) - No MC labels available, please check!!!\n"); |
257 | return kMCBadLabel; | |
b5ef1905 | 258 | } |
259 | ||
3e57baab | 260 | TObjArray* arrayCluster = 0; |
261 | if ( calorimeter == "EMCAL" ) arrayCluster = reader->GetEMCALClusters(); | |
262 | else if ( calorimeter == "PHOS" ) arrayCluster= reader->GetPHOSClusters(); | |
263 | ||
3c75ddf7 | 264 | Int_t labels[]={label}; |
265 | ||
266 | //Select where the information is, ESD-galice stack or AOD mcparticles branch | |
267 | if(reader->ReadStack()){ | |
3e57baab | 268 | tag = CheckOriginInStack(labels, 1,reader->GetStack(),arrayCluster); |
3c75ddf7 | 269 | } |
270 | else if(reader->ReadAODMCParticles()){ | |
3e57baab | 271 | tag = CheckOriginInAOD(labels, 1,reader->GetAODMCParticles(),arrayCluster); |
3c75ddf7 | 272 | } |
273 | ||
274 | return tag ; | |
7cd4e982 | 275 | } |
276 | ||
3e57baab | 277 | //__________________________________________________________________________________________ |
278 | Int_t AliMCAnalysisUtils::CheckOriginInStack(const Int_t *labels, Int_t nlabels, | |
279 | AliStack* stack, const TObjArray* arrayCluster) | |
840124fe | 280 | { |
7cd4e982 | 281 | // Play with the MC stack if available. Tag particles depending on their origin. |
282 | // Do same things as in CheckOriginInAOD but different input. | |
556e55b0 | 283 | |
284 | //generally speaking, label is the MC label of a reconstructed | |
285 | //entity (track, cluster, etc) for which we want to know something | |
286 | //about its heritage, but one can also use it directly with stack | |
287 | //particles not connected to reconstructed entities | |
f8006433 | 288 | |
36769d30 | 289 | if(!stack) |
290 | { | |
1cd71065 | 291 | if (fDebug >=0) |
f8006433 | 292 | printf("AliMCAnalysisUtils::CheckOriginInStack() - Stack is not available, check analysis settings in configuration file, STOP!!\n"); |
293 | return -1; | |
7cd4e982 | 294 | } |
f8006433 | 295 | |
7cd4e982 | 296 | Int_t tag = 0; |
902aa95c | 297 | Int_t label=labels[0];//Most significant particle contributing to the cluster |
3c75ddf7 | 298 | |
36769d30 | 299 | if(label >= 0 && label < stack->GetNtrack()) |
300 | { | |
556e55b0 | 301 | //MC particle of interest is the "mom" of the entity |
7cd4e982 | 302 | TParticle * mom = stack->Particle(label); |
eb3e1b50 | 303 | Int_t iMom = label; |
304 | Int_t mPdgSign = mom->GetPdgCode(); | |
305 | Int_t mPdg = TMath::Abs(mPdgSign); | |
306 | Int_t mStatus = mom->GetStatusCode() ; | |
307 | Int_t iParent = mom->GetFirstMother() ; | |
36769d30 | 308 | |
0d360687 | 309 | if(fDebug > 0 && label < 8 && fMCGenerator != kBoxLike) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent); |
7cd4e982 | 310 | |
556e55b0 | 311 | //GrandParent of the entity |
d7c10d78 | 312 | TParticle * parent = NULL; |
7cd4e982 | 313 | Int_t pPdg = -1; |
314 | Int_t pStatus =-1; | |
36769d30 | 315 | if(iParent >= 0) |
316 | { | |
7cd4e982 | 317 | parent = stack->Particle(iParent); |
36769d30 | 318 | |
319 | if(parent) | |
320 | { | |
95fb4ab7 | 321 | pPdg = TMath::Abs(parent->GetPdgCode()); |
322 | pStatus = parent->GetStatusCode(); | |
323 | } | |
7cd4e982 | 324 | } |
325 | else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Parent with label %d\n",iParent); | |
8dacfd76 | 326 | |
36769d30 | 327 | if(fDebug > 2 ) |
328 | { | |
3c75ddf7 | 329 | printf("AliMCAnalysisUtils::CheckOriginInStack() - Cluster most contributing mother and its parent: \n"); |
330 | printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus); | |
331 | printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus); | |
f8006433 | 332 | } |
902aa95c | 333 | |
556e55b0 | 334 | //Check if "mother" of entity is converted, if not, get the first non converted mother |
36769d30 | 335 | if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && mStatus == 0) |
336 | { | |
8dacfd76 | 337 | SetTagBit(tag,kMCConversion); |
36769d30 | 338 | |
8dacfd76 | 339 | //Check if the mother is photon or electron with status not stable |
36769d30 | 340 | while ((pPdg == 22 || pPdg == 11) && mStatus != 1) |
341 | { | |
f8006433 | 342 | //Mother |
eb3e1b50 | 343 | iMom = mom->GetFirstMother(); |
344 | mom = stack->Particle(iMom); | |
345 | mPdgSign = mom->GetPdgCode(); | |
346 | mPdg = TMath::Abs(mPdgSign); | |
347 | mStatus = mom->GetStatusCode() ; | |
348 | iParent = mom->GetFirstMother() ; | |
36769d30 | 349 | |
f8006433 | 350 | if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent); |
351 | ||
352 | //GrandParent | |
36769d30 | 353 | if(iParent >= 0) |
354 | { | |
f8006433 | 355 | parent = stack->Particle(iParent); |
36769d30 | 356 | if(parent) |
357 | { | |
95fb4ab7 | 358 | pPdg = TMath::Abs(parent->GetPdgCode()); |
359 | pStatus = parent->GetStatusCode(); | |
360 | } | |
f8006433 | 361 | } |
0e453907 | 362 | else {// in case of gun/box simulations |
363 | pPdg = 0; | |
364 | pStatus = 0; | |
365 | break; | |
366 | } | |
36769d30 | 367 | }//while |
368 | ||
369 | if(fDebug > 2 ) | |
370 | { | |
f8006433 | 371 | printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted photon/electron: \n"); |
372 | printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus); | |
373 | printf("\t Parent label %d, pdg %d, status %d\n",iParent, pPdg, pStatus); | |
374 | } | |
375 | ||
7cd4e982 | 376 | }//mother and parent are electron or photon and have status 0 |
36769d30 | 377 | else if((mPdg == 22 || mPdg == 11) && mStatus == 0) |
378 | { | |
902aa95c | 379 | //Still a conversion but only one electron/photon generated. Just from hadrons but not decays. |
380 | if(pPdg == 2112 || pPdg == 211 || pPdg == 321 || | |
36769d30 | 381 | pPdg == 2212 || pPdg == 130 || pPdg == 13 ) |
382 | { | |
f8006433 | 383 | SetTagBit(tag,kMCConversion); |
eb3e1b50 | 384 | iMom = mom->GetFirstMother(); |
385 | mom = stack->Particle(iMom); | |
386 | mPdgSign = mom->GetPdgCode(); | |
387 | mPdg = TMath::Abs(mPdgSign); | |
f8006433 | 388 | |
36769d30 | 389 | if(fDebug > 2 ) |
390 | { | |
f8006433 | 391 | printf("AliMCAnalysisUtils::CheckOriginInStack() - Converted hadron: \n"); |
392 | printf("\t Mother label %d, pdg %d, status %d\n",iMom, mPdg, mStatus); | |
393 | } | |
394 | }//hadron converted | |
395 | ||
8dacfd76 | 396 | //Comment for the next lines, we do not check the parent of the hadron for the moment. |
397 | //iParent = mom->GetFirstMother() ; | |
398 | //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - Mother is parton %d\n",iParent); | |
399 | ||
400 | //GrandParent | |
401 | //if(iParent >= 0){ | |
402 | // parent = stack->Particle(iParent); | |
403 | // pPdg = TMath::Abs(parent->GetPdgCode()); | |
404 | //} | |
405 | } | |
902aa95c | 406 | // conversion into electrons/photons checked |
8dacfd76 | 407 | |
7cd4e982 | 408 | //first check for typical charged particles |
eb3e1b50 | 409 | if (mPdg == 13) SetTagBit(tag,kMCMuon); |
410 | else if(mPdg == 211) SetTagBit(tag,kMCPion); | |
411 | else if(mPdg == 321) SetTagBit(tag,kMCKaon); | |
412 | else if(mPdgSign == 2212) SetTagBit(tag,kMCProton); | |
413 | else if(mPdgSign == -2212) SetTagBit(tag,kMCAntiProton); | |
414 | else if(mPdgSign == 2112) SetTagBit(tag,kMCNeutron); | |
415 | else if(mPdgSign == -2112) SetTagBit(tag,kMCAntiNeutron); | |
840124fe | 416 | |
8dacfd76 | 417 | //check for pi0 and eta (shouldn't happen unless their decays were turned off) |
36769d30 | 418 | else if(mPdg == 111) |
419 | { | |
4d52b003 | 420 | |
f8006433 | 421 | SetTagBit(tag,kMCPi0Decay); |
36769d30 | 422 | |
4d52b003 | 423 | if(fDebug > 2 ) |
424 | printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly pi0, not decayed by generator \n"); | |
36769d30 | 425 | |
f8006433 | 426 | CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCPi0 if 2 gammas in same cluster |
3e57baab | 427 | |
f8006433 | 428 | } |
36769d30 | 429 | else if(mPdg == 221) |
430 | { | |
f8006433 | 431 | SetTagBit(tag,kMCEtaDecay); |
36769d30 | 432 | |
4d52b003 | 433 | if(fDebug > 2 ) |
434 | printf("AliMCAnalysisUtils::CheckOriginInStack() - First mother is directly eta, not decayed by generator \n"); | |
36769d30 | 435 | |
f8006433 | 436 | CheckOverlapped2GammaDecay(labels,nlabels, iMom, stack, tag); //set to kMCEta if 2 gammas in same cluster |
437 | } | |
8dacfd76 | 438 | //Photons |
36769d30 | 439 | else if(mPdg == 22) |
440 | { | |
7cd4e982 | 441 | SetTagBit(tag,kMCPhoton); |
4d52b003 | 442 | |
443 | if(pPdg == 111) | |
444 | { | |
445 | SetTagBit(tag,kMCPi0Decay); | |
446 | ||
447 | if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA pi0 decay photon, parent pi0 with status 11 \n"); | |
448 | ||
449 | CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag); //set to kMCPi0 if 2 gammas in same cluster | |
3e57baab | 450 | // In case it did not merge, check if the decay companion is lost |
451 | if(!CheckTagBit(tag, kMCPi0) && !CheckTagBit(tag,kMCDecayPairInCalo)) | |
452 | CheckLostDecayPair(arrayCluster,iMom, iParent, stack, tag); | |
453 | ||
454 | //printf("Bit set is Merged %d, Pair in calo %d, Lost %d\n",CheckTagBit(tag, kMCPi0),CheckTagBit(tag,kMCDecayPairInCalo),CheckTagBit(tag,kMCDecayPairLost)); | |
4d52b003 | 455 | } |
456 | else if (pPdg == 221) | |
457 | { | |
458 | SetTagBit(tag, kMCEtaDecay); | |
459 | ||
460 | if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInStack() - PYTHIA eta decay photon, parent pi0 with status 11 \n"); | |
461 | ||
462 | CheckOverlapped2GammaDecay(labels,nlabels, iParent, stack, tag);//set to kMCEta if 2 gammas in same cluster | |
3e57baab | 463 | // In case it did not merge, check if the decay companion is lost |
464 | if(!CheckTagBit(tag, kMCEta) && !CheckTagBit(tag,kMCDecayPairInCalo)) | |
465 | CheckLostDecayPair(arrayCluster,iMom, iParent, stack, tag); | |
4d52b003 | 466 | } |
467 | else if(mStatus == 1) | |
36769d30 | 468 | { //undecayed particle |
0d360687 | 469 | if(fMCGenerator == kPythia) |
36769d30 | 470 | { |
471 | if(iParent < 8 && iParent > 5) | |
472 | {//outgoing partons | |
f8006433 | 473 | if(pPdg == 22) SetTagBit(tag,kMCPrompt); |
36769d30 | 474 | else SetTagBit(tag,kMCFragmentation); |
f8006433 | 475 | }//Outgoing partons |
36769d30 | 476 | else if(iParent <= 5) |
477 | { | |
f8006433 | 478 | SetTagBit(tag, kMCISR); //Initial state radiation |
479 | } | |
4d52b003 | 480 | else SetTagBit(tag,kMCUnknown); |
481 | }//PYTHIA | |
482 | ||
0d360687 | 483 | else if(fMCGenerator == kHerwig) |
36769d30 | 484 | { |
485 | if(pStatus < 197) | |
486 | {//Not decay | |
487 | while(1) | |
488 | { | |
489 | if(parent) | |
490 | { | |
95fb4ab7 | 491 | if(parent->GetFirstMother()<=5) break; |
492 | iParent = parent->GetFirstMother(); | |
493 | parent=stack->Particle(iParent); | |
494 | pStatus= parent->GetStatusCode(); | |
495 | pPdg = TMath::Abs(parent->GetPdgCode()); | |
496 | } else break; | |
f8006433 | 497 | }//Look for the parton |
498 | ||
36769d30 | 499 | if(iParent < 8 && iParent > 5) |
500 | { | |
f8006433 | 501 | if(pPdg == 22) SetTagBit(tag,kMCPrompt); |
36769d30 | 502 | else SetTagBit(tag,kMCFragmentation); |
f8006433 | 503 | } |
504 | else SetTagBit(tag,kMCISR);//Initial state radiation | |
505 | }//Not decay | |
4d52b003 | 506 | else SetTagBit(tag,kMCUnknown); |
f8006433 | 507 | }//HERWIG |
4d52b003 | 508 | } |
509 | else SetTagBit(tag,kMCOtherDecay); | |
510 | ||
7cd4e982 | 511 | }//Mother Photon |
8dacfd76 | 512 | |
7cd4e982 | 513 | //Electron check. Where did that electron come from? |
514 | else if(mPdg == 11){ //electron | |
36769d30 | 515 | if(pPdg == 11 && parent) |
516 | { | |
e495fc0d | 517 | Int_t iGrandma = parent->GetFirstMother(); |
36769d30 | 518 | if(iGrandma >= 0) |
519 | { | |
e495fc0d | 520 | TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother |
521 | Int_t gPdg = TMath::Abs(gma->GetPdgCode()); | |
f8006433 | 522 | |
36769d30 | 523 | if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson |
e495fc0d | 524 | else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson |
525 | } | |
526 | } | |
36769d30 | 527 | |
528 | SetTagBit(tag,kMCElectron); | |
529 | ||
902aa95c | 530 | if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInStack() - Checking ancestors of electrons\n"); |
36769d30 | 531 | |
3e57baab | 532 | if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); SetTagBit(tag, kMCDecayDalitz) ; } //Pi0 Dalitz decay |
533 | else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); SetTagBit(tag, kMCDecayDalitz) ; } //Eta Dalitz decay | |
d409d6b9 | 534 | else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB); } //b-->e decay |
b3fdfed5 | 535 | else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000)) { //check charm decay |
36769d30 | 536 | if(parent) |
537 | { | |
b3fdfed5 | 538 | Int_t iGrandma = parent->GetFirstMother(); |
36769d30 | 539 | if(iGrandma >= 0) |
540 | { | |
b3fdfed5 | 541 | TParticle* gma = (TParticle*)stack->Particle(iGrandma); //get mother of charm |
542 | Int_t gPdg = TMath::Abs(gma->GetPdgCode()); | |
543 | if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e | |
544 | else SetTagBit(tag,kMCEFromC); //c-->e | |
36769d30 | 545 | } |
546 | else SetTagBit(tag,kMCEFromC); //c-->e | |
b3fdfed5 | 547 | }//parent |
36769d30 | 548 | } |
549 | else | |
550 | { | |
f8006433 | 551 | //if it is not from any of the above, where is it from? |
552 | if(pPdg > 10000) SetTagBit(tag,kMCUnknown); | |
36769d30 | 553 | |
f8006433 | 554 | else SetTagBit(tag,kMCOtherDecay); |
36769d30 | 555 | |
95fb4ab7 | 556 | if(fDebug > 0 && parent) printf("AliMCAnalysisUtils::CheckOriginInStack() - Status %d Electron from other origin: %s (pPdg = %d) %s (mpdg = %d)\n",mStatus,parent->GetName(),pPdg,mom->GetName(),mPdg); |
d409d6b9 | 557 | } |
7cd4e982 | 558 | }//electron check |
559 | //Cluster was made by something else | |
36769d30 | 560 | else |
561 | { | |
562 | if(fDebug > 0) | |
563 | printf("AliMCAnalysisUtils::CheckOriginInStack() - \tSetting kMCUnknown for cluster from %s (pdg = %d, Parent pdg = %d)\n", | |
564 | mom->GetName(),mPdg,pPdg); | |
565 | ||
7cd4e982 | 566 | SetTagBit(tag,kMCUnknown); |
567 | } | |
568 | }//Good label value | |
36769d30 | 569 | else |
570 | {// Bad label | |
902aa95c | 571 | if(label < 0 && (fDebug >= 0)) |
f8006433 | 572 | printf("AliMCAnalysisUtils::CheckOriginInStack() *** bad label or no stack ***: label %d \n", label); |
36769d30 | 573 | |
574 | if(label >= stack->GetNtrack() && (fDebug >= 0)) | |
f8006433 | 575 | printf("AliMCAnalysisUtils::CheckOriginInStack() *** large label ***: label %d, n tracks %d \n", label, stack->GetNtrack()); |
36769d30 | 576 | |
7cd4e982 | 577 | SetTagBit(tag,kMCUnknown); |
578 | }//Bad label | |
579 | ||
580 | return tag; | |
581 | ||
582 | } | |
583 | ||
584 | ||
3e57baab | 585 | //________________________________________________________________________________________________________ |
8a2dbbff | 586 | Int_t AliMCAnalysisUtils::CheckOriginInAOD(const Int_t *labels, Int_t nlabels, |
3e57baab | 587 | const TClonesArray *mcparticles, const TObjArray* arrayCluster) |
840124fe | 588 | { |
7cd4e982 | 589 | // Play with the MCParticles in AOD if available. Tag particles depending on their origin. |
590 | // Do same things as in CheckOriginInStack but different input. | |
0d360687 | 591 | |
36769d30 | 592 | if(!mcparticles) |
593 | { | |
1cd71065 | 594 | if(fDebug >= 0) |
f8006433 | 595 | printf("AliMCAnalysisUtils::CheckOriginInAOD() - AODMCParticles is not available, check analysis settings in configuration file!!\n"); |
596 | return -1; | |
7cd4e982 | 597 | } |
7cd4e982 | 598 | |
599 | Int_t tag = 0; | |
902aa95c | 600 | Int_t label=labels[0];//Most significant particle contributing to the cluster |
f8006433 | 601 | |
7cd4e982 | 602 | Int_t nprimaries = mcparticles->GetEntriesFast(); |
36769d30 | 603 | if(label >= 0 && label < nprimaries) |
604 | { | |
7cd4e982 | 605 | //Mother |
606 | AliAODMCParticle * mom = (AliAODMCParticle *) mcparticles->At(label); | |
eb3e1b50 | 607 | Int_t iMom = label; |
608 | Int_t mPdgSign = mom->GetPdgCode(); | |
609 | Int_t mPdg = TMath::Abs(mPdgSign); | |
610 | Int_t iParent = mom->GetMother() ; | |
36769d30 | 611 | |
0d360687 | 612 | if(fDebug > 0 && label < 8 && fMCGenerator != kBoxLike) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent); |
7cd4e982 | 613 | |
614 | //GrandParent | |
95fb4ab7 | 615 | AliAODMCParticle * parent = NULL ; |
7cd4e982 | 616 | Int_t pPdg = -1; |
36769d30 | 617 | if(iParent >= 0) |
618 | { | |
7cd4e982 | 619 | parent = (AliAODMCParticle *) mcparticles->At(iParent); |
620 | pPdg = TMath::Abs(parent->GetPdgCode()); | |
621 | } | |
622 | else if(fDebug > 0 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Parent with label %d\n",iParent); | |
f8006433 | 623 | |
36769d30 | 624 | if(fDebug > 2 ) |
625 | { | |
3c75ddf7 | 626 | printf("AliMCAnalysisUtils::CheckOriginInAOD() - Cluster most contributing mother and its parent: \n"); |
36769d30 | 627 | |
628 | printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n", | |
629 | iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary()); | |
630 | ||
3c75ddf7 | 631 | if(parent) |
36769d30 | 632 | printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n", |
633 | iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary()); | |
f8006433 | 634 | } |
902aa95c | 635 | |
8dacfd76 | 636 | //Check if mother is converted, if not, get the first non converted mother |
36769d30 | 637 | if((mPdg == 22 || mPdg == 11) && (pPdg == 22 || pPdg == 11) && !mom->IsPrimary()) |
638 | { | |
8dacfd76 | 639 | SetTagBit(tag,kMCConversion); |
36769d30 | 640 | |
8dacfd76 | 641 | //Check if the mother is photon or electron with status not stable |
36769d30 | 642 | while ((pPdg == 22 || pPdg == 11) && !mom->IsPhysicalPrimary()) |
643 | { | |
f8006433 | 644 | //Mother |
eb3e1b50 | 645 | iMom = mom->GetMother(); |
646 | mom = (AliAODMCParticle *) mcparticles->At(iMom); | |
647 | mPdgSign = mom->GetPdgCode(); | |
648 | mPdg = TMath::Abs(mPdgSign); | |
649 | iParent = mom->GetMother() ; | |
f8006433 | 650 | if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent); |
651 | ||
652 | //GrandParent | |
36769d30 | 653 | if(iParent >= 0 && parent) |
654 | { | |
f8006433 | 655 | parent = (AliAODMCParticle *) mcparticles->At(iParent); |
656 | pPdg = TMath::Abs(parent->GetPdgCode()); | |
657 | } | |
658 | // printf("\t While Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary()); | |
659 | // printf("\t While Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary()); | |
660 | ||
902aa95c | 661 | }//while |
f8006433 | 662 | |
36769d30 | 663 | if(fDebug > 2 ) |
664 | { | |
f8006433 | 665 | printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted photon/electron : \n"); |
36769d30 | 666 | |
667 | printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n", | |
668 | iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary()); | |
669 | ||
95fb4ab7 | 670 | if(parent) |
36769d30 | 671 | printf("\t Parent label %d, pdg %d, Primary? %d, Physical Primary? %d\n", |
672 | iParent, pPdg, parent->IsPrimary(), parent->IsPhysicalPrimary()); | |
f8006433 | 673 | } |
674 | ||
8dacfd76 | 675 | }//mother and parent are electron or photon and have status 0 and parent is photon or electron |
36769d30 | 676 | else if((mPdg == 22 || mPdg == 11) && !mom->IsPrimary()) |
677 | { | |
8dacfd76 | 678 | //Still a conversion but only one electron/photon generated. Just from hadrons |
902aa95c | 679 | if(pPdg == 2112 || pPdg == 211 || pPdg == 321 || |
36769d30 | 680 | pPdg == 2212 || pPdg == 130 || pPdg == 13 ) |
681 | { | |
f8006433 | 682 | SetTagBit(tag,kMCConversion); |
eb3e1b50 | 683 | iMom = mom->GetMother(); |
684 | mom = (AliAODMCParticle *) mcparticles->At(iMom); | |
685 | mPdgSign = mom->GetPdgCode(); | |
686 | mPdg = TMath::Abs(mPdgSign); | |
f8006433 | 687 | |
36769d30 | 688 | if(fDebug > 2 ) |
689 | { | |
f8006433 | 690 | printf("AliMCAnalysisUtils::CheckOriginInAOD() - Converted hadron : \n"); |
691 | printf("\t Mother label %d, pdg %d, Primary? %d, Physical Primary? %d\n",iMom, mPdg, mom->IsPrimary(), mom->IsPhysicalPrimary()); | |
692 | } | |
693 | }//hadron converted | |
694 | ||
8dacfd76 | 695 | //Comment for next lines, we do not check the parent of the hadron for the moment. |
696 | //iParent = mom->GetMother() ; | |
697 | //if(fDebug > 0 && label < 8 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Mother is parton %d\n",iParent); | |
698 | ||
699 | //GrandParent | |
700 | //if(iParent >= 0){ | |
701 | // parent = (AliAODMCParticle *) mcparticles->At(iParent); | |
702 | // pPdg = TMath::Abs(parent->GetPdgCode()); | |
703 | //} | |
704 | } | |
705 | ||
3c75ddf7 | 706 | //printf("Final mother mPDG %d\n",mPdg); |
c5693f62 | 707 | |
8dacfd76 | 708 | // conversion into electrons/photons checked |
709 | ||
710 | //first check for typical charged particles | |
eb3e1b50 | 711 | if (mPdg == 13) SetTagBit(tag,kMCMuon); |
712 | else if(mPdg == 211) SetTagBit(tag,kMCPion); | |
713 | else if(mPdg == 321) SetTagBit(tag,kMCKaon); | |
714 | else if(mPdgSign == 2212) SetTagBit(tag,kMCProton); | |
715 | else if(mPdgSign == 2112) SetTagBit(tag,kMCNeutron); | |
716 | else if(mPdgSign == -2212) SetTagBit(tag,kMCAntiProton); | |
717 | else if(mPdgSign == -2112) SetTagBit(tag,kMCAntiNeutron); | |
718 | ||
8dacfd76 | 719 | //check for pi0 and eta (shouldn't happen unless their decays were turned off) |
36769d30 | 720 | else if(mPdg == 111) |
721 | { | |
f8006433 | 722 | SetTagBit(tag,kMCPi0Decay); |
36769d30 | 723 | |
f8006433 | 724 | if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly pi0, not decayed by generator \n"); |
36769d30 | 725 | |
f8006433 | 726 | CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster |
727 | } | |
36769d30 | 728 | else if(mPdg == 221) |
729 | { | |
f8006433 | 730 | SetTagBit(tag,kMCEtaDecay); |
36769d30 | 731 | |
f8006433 | 732 | if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - First mother is directly eta, not decayed by generator \n"); |
36769d30 | 733 | |
f8006433 | 734 | CheckOverlapped2GammaDecay(labels,nlabels, iMom, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster |
735 | } | |
8dacfd76 | 736 | //Photons |
36769d30 | 737 | else if(mPdg == 22) |
738 | { | |
7cd4e982 | 739 | SetTagBit(tag,kMCPhoton); |
4d52b003 | 740 | |
741 | if(pPdg == 111) | |
742 | { | |
743 | SetTagBit(tag,kMCPi0Decay); | |
744 | ||
745 | if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator pi0 decay photon \n"); | |
746 | ||
747 | CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCPi0 if 2 gammas in same cluster | |
3e57baab | 748 | // In case it did not merge, check if the decay companion is lost |
749 | if(!CheckTagBit(tag, kMCPi0) && !CheckTagBit(tag,kMCDecayPairInCalo) && !CheckTagBit(tag,kMCDecayPairLost)) | |
750 | CheckLostDecayPair(arrayCluster,iMom, iParent, mcparticles, tag); | |
751 | ||
752 | //printf("Bit set is Merged %d, Pair in calo %d, Lost %d\n",CheckTagBit(tag, kMCPi0),CheckTagBit(tag,kMCDecayPairInCalo),CheckTagBit(tag,kMCDecayPairLost)); | |
4d52b003 | 753 | } |
754 | else if (pPdg == 221) | |
755 | { | |
756 | SetTagBit(tag, kMCEtaDecay); | |
757 | ||
758 | if(fDebug > 2 ) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Generator eta decay photon \n"); | |
759 | ||
760 | CheckOverlapped2GammaDecay(labels,nlabels, iParent, mcparticles, tag); //set to kMCEta if 2 gammas in same cluster | |
3e57baab | 761 | // In case it did not merge, check if the decay companion is lost |
762 | if(!CheckTagBit(tag, kMCEta) && !CheckTagBit(tag,kMCDecayPairInCalo)) | |
763 | CheckLostDecayPair(arrayCluster,iMom, iParent, mcparticles, tag); | |
4d52b003 | 764 | } |
0d360687 | 765 | else if( mom->IsPhysicalPrimary() && ( fMCGenerator == kPythia || fMCGenerator == kHerwig ) ) //undecayed particle |
3c75ddf7 | 766 | { |
36769d30 | 767 | if(iParent < 8 && iParent > 5 ) |
768 | {//outgoing partons | |
3c75ddf7 | 769 | if(pPdg == 22) SetTagBit(tag,kMCPrompt); |
770 | else SetTagBit(tag,kMCFragmentation); | |
771 | }//Outgoing partons | |
0d360687 | 772 | else if( iParent <= 5 && ( fMCGenerator == kPythia || fMCGenerator == kHerwig ) ) |
36769d30 | 773 | { |
3c75ddf7 | 774 | SetTagBit(tag, kMCISR); //Initial state radiation |
775 | } | |
4d52b003 | 776 | else SetTagBit(tag,kMCUnknown); |
840124fe | 777 | }//Physical primary |
4d52b003 | 778 | else SetTagBit(tag,kMCOtherDecay); |
779 | ||
7cd4e982 | 780 | }//Mother Photon |
781 | ||
782 | //Electron check. Where did that electron come from? | |
36769d30 | 783 | else if(mPdg == 11) |
784 | { //electron | |
785 | if(pPdg == 11 && parent) | |
786 | { | |
e495fc0d | 787 | Int_t iGrandma = parent->GetMother(); |
36769d30 | 788 | if(iGrandma >= 0) |
789 | { | |
e495fc0d | 790 | AliAODMCParticle* gma = (AliAODMCParticle*)mcparticles->At(iGrandma); |
791 | Int_t gPdg = TMath::Abs(gma->GetPdgCode()); | |
f8006433 | 792 | |
36769d30 | 793 | if (gPdg == 23) { SetTagBit(tag,kMCZDecay); } //parent is Z-boson |
e495fc0d | 794 | else if (gPdg == 24) { SetTagBit(tag,kMCWDecay); } //parent is W-boson |
795 | } | |
796 | } | |
36769d30 | 797 | |
798 | SetTagBit(tag,kMCElectron); | |
799 | ||
3e57baab | 800 | if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Checking ancestors of electrons"); |
801 | ||
802 | if (pPdg == 111) { SetTagBit(tag,kMCPi0Decay); SetTagBit(tag,kMCDecayDalitz);} //Pi0 Dalitz decay | |
803 | else if (pPdg == 221) { SetTagBit(tag,kMCEtaDecay); SetTagBit(tag,kMCDecayDalitz);} //Eta Dalitz decay | |
d409d6b9 | 804 | else if((499 < pPdg && pPdg < 600)||(4999 < pPdg && pPdg < 6000)) { SetTagBit(tag,kMCEFromB);} //b-hadron decay |
36769d30 | 805 | else if((399 < pPdg && pPdg < 500)||(3999 < pPdg && pPdg < 5000)) |
806 | { //c-hadron decay check | |
807 | if(parent) | |
808 | { | |
b3fdfed5 | 809 | Int_t iGrandma = parent->GetMother(); |
36769d30 | 810 | if(iGrandma >= 0) |
811 | { | |
b3fdfed5 | 812 | AliAODMCParticle* gma = (AliAODMCParticle*)mcparticles->At(iGrandma); //charm's mother |
813 | Int_t gPdg = TMath::Abs(gma->GetPdgCode()); | |
814 | if((499 < gPdg && gPdg < 600)||(4999 < gPdg && gPdg < 6000)) SetTagBit(tag,kMCEFromCFromB); //b-->c-->e decay | |
815 | else SetTagBit(tag,kMCEFromC); //c-hadron decay | |
36769d30 | 816 | } |
817 | else SetTagBit(tag,kMCEFromC); //c-hadron decay | |
b3fdfed5 | 818 | }//parent |
36769d30 | 819 | } else |
820 | { //prompt or other decay | |
f8006433 | 821 | TParticlePDG* foo = TDatabasePDG::Instance()->GetParticle(pPdg); |
822 | TParticlePDG* foo1 = TDatabasePDG::Instance()->GetParticle(mPdg); | |
36769d30 | 823 | |
824 | if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - Electron from other origin: %s (pPdg = %d) %s (mPdg = %d)\n", | |
825 | foo->GetName(), pPdg,foo1->GetName(),mPdg); | |
826 | ||
f8006433 | 827 | if(pPdg > 10000) SetTagBit(tag,kMCUnknown); |
36769d30 | 828 | else SetTagBit(tag,kMCOtherDecay); |
d409d6b9 | 829 | } |
8dacfd76 | 830 | }//electron check |
7cd4e982 | 831 | //cluster was made by something else |
36769d30 | 832 | else |
833 | { | |
834 | if(fDebug > 0) printf("AliMCAnalysisUtils::CheckOriginInAOD() - \tSetting kMCUnknown for cluster with pdg = %d, Parent pdg = %d\n", | |
835 | mPdg,pPdg); | |
7cd4e982 | 836 | SetTagBit(tag,kMCUnknown); |
837 | } | |
838 | }//Good label value | |
36769d30 | 839 | else |
840 | {//Bad label | |
902aa95c | 841 | if(label < 0 && (fDebug >= 0) ) |
f8006433 | 842 | printf("AliMCAnalysisUtils::CheckOriginInAOD() *** bad label or no mcparticles ***: label %d \n", label); |
36769d30 | 843 | |
844 | if(label >= mcparticles->GetEntriesFast() && (fDebug >= 0) ) | |
f8006433 | 845 | printf("AliMCAnalysisUtils::CheckOriginInAOD() *** large label ***: label %d, n tracks %d \n", label, mcparticles->GetEntriesFast()); |
36769d30 | 846 | |
7cd4e982 | 847 | SetTagBit(tag,kMCUnknown); |
f8006433 | 848 | |
7cd4e982 | 849 | }//Bad label |
850 | ||
851 | return tag; | |
852 | ||
853 | } | |
854 | ||
8a2dbbff | 855 | //_________________________________________________________________________________________ |
856 | void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels, Int_t nlabels, | |
857 | Int_t mesonIndex, AliStack *stack, | |
c5693f62 | 858 | Int_t &tag) |
840124fe | 859 | { |
0d360687 | 860 | // Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in stack. |
3c75ddf7 | 861 | |
862 | if(labels[0] < 0 || labels[0] > stack->GetNtrack() || nlabels <= 1) { | |
36769d30 | 863 | if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Exit : label[0] %d, n primaries %d, nlabels %d \n", |
f8006433 | 864 | labels[0],stack->GetNtrack(), nlabels); |
3c75ddf7 | 865 | return; |
866 | } | |
867 | ||
868 | TParticle * meson = stack->Particle(mesonIndex); | |
869 | Int_t mesonPdg = meson->GetPdgCode(); | |
36769d30 | 870 | if(mesonPdg!=111 && mesonPdg!=221) |
871 | { | |
872 | printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Wrong pi0/eta PDG : %d \n",mesonPdg); | |
3c75ddf7 | 873 | return; |
874 | } | |
875 | ||
36769d30 | 876 | if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - %s, label %d\n",meson->GetName(), mesonIndex); |
3c75ddf7 | 877 | |
878 | //Check if meson decayed into 2 daughters or if both were kept. | |
36769d30 | 879 | if(meson->GetNDaughters() != 2) |
880 | { | |
3c75ddf7 | 881 | if(fDebug > 2) |
36769d30 | 882 | printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters()); |
3c75ddf7 | 883 | return; |
884 | } | |
885 | ||
886 | //Get the daughters | |
887 | Int_t iPhoton0 = meson->GetDaughter(0); | |
888 | Int_t iPhoton1 = meson->GetDaughter(1); | |
889 | TParticle *photon0 = stack->Particle(iPhoton0); | |
890 | TParticle *photon1 = stack->Particle(iPhoton1); | |
891 | ||
892 | //Check if both daughters are photons | |
36769d30 | 893 | if(photon0->GetPdgCode() != 22 || photon1->GetPdgCode()!=22) |
894 | { | |
3c75ddf7 | 895 | if(fDebug > 2) |
36769d30 | 896 | printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Not overalapped. PDG: daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode()); |
3c75ddf7 | 897 | return; |
898 | } | |
899 | ||
900 | if(fDebug > 2) | |
36769d30 | 901 | printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1); |
3c75ddf7 | 902 | |
903 | //Check if both photons contribute to the cluster | |
904 | Bool_t okPhoton0 = kFALSE; | |
905 | Bool_t okPhoton1 = kFALSE; | |
906 | ||
36769d30 | 907 | if(fDebug > 3) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Labels loop:\n"); |
908 | ||
909 | Bool_t conversion = kFALSE; | |
3c75ddf7 | 910 | |
36769d30 | 911 | for(Int_t i = 0; i < nlabels; i++) |
912 | { | |
3c75ddf7 | 913 | if(fDebug > 3) printf("\t at begin:label %d/%d: %d, ok? photon1 %d, photon2 %d\n", i+1, nlabels, labels[i], okPhoton0, okPhoton1); |
914 | ||
915 | //If we already found both, break the loop | |
916 | if(okPhoton0 && okPhoton1) break; | |
917 | ||
36769d30 | 918 | Int_t index = labels[i]; |
919 | if (iPhoton0 == index) | |
920 | { | |
3c75ddf7 | 921 | okPhoton0 = kTRUE; |
922 | continue; | |
923 | } | |
36769d30 | 924 | else if (iPhoton1 == index) |
925 | { | |
3c75ddf7 | 926 | okPhoton1 = kTRUE; |
927 | continue; | |
928 | } | |
929 | ||
930 | //Trace back the mother in case it was a conversion | |
931 | ||
36769d30 | 932 | if(index >= stack->GetNtrack()) |
933 | { | |
934 | printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) Particle index %d larger than size of list %d !!\n", | |
935 | index,stack->GetNtrack()); | |
3c75ddf7 | 936 | continue; |
937 | } | |
938 | ||
939 | TParticle * daught = stack->Particle(index); | |
940 | Int_t tmpindex = daught->GetFirstMother(); | |
941 | if(fDebug > 3) printf("\t Conversion? : mother %d\n",tmpindex); | |
36769d30 | 942 | while(tmpindex>=0) |
943 | { | |
3c75ddf7 | 944 | //MC particle of interest is the mother |
945 | if(fDebug > 3) printf("\t \t parent index %d\n",tmpindex); | |
946 | daught = stack->Particle(tmpindex); | |
36769d30 | 947 | if (iPhoton0 == tmpindex) |
948 | { | |
949 | conversion = kTRUE; | |
950 | okPhoton0 = kTRUE; | |
c5693f62 | 951 | break; |
3c75ddf7 | 952 | } |
36769d30 | 953 | else if (iPhoton1 == tmpindex) |
954 | { | |
955 | conversion = kTRUE; | |
956 | okPhoton1 = kTRUE; | |
c5693f62 | 957 | break; |
3c75ddf7 | 958 | } |
36769d30 | 959 | |
3c75ddf7 | 960 | tmpindex = daught->GetFirstMother(); |
36769d30 | 961 | |
3c75ddf7 | 962 | }//While to check if pi0/eta daughter was one of these contributors to the cluster |
963 | ||
964 | if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0) | |
36769d30 | 965 | printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Something happens, first label should be from a photon decay!\n"); |
3c75ddf7 | 966 | |
967 | }//loop on list of labels | |
968 | ||
969 | //If both photons contribute tag as the corresponding meson. | |
36769d30 | 970 | if(okPhoton0 && okPhoton1) |
971 | { | |
3c75ddf7 | 972 | if(fDebug > 2) |
36769d30 | 973 | printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - %s OVERLAPPED DECAY \n", meson->GetName()); |
974 | ||
975 | if(!CheckTagBit(tag,kMCConversion) && conversion) SetTagBit(tag,kMCConversion) ; | |
3c75ddf7 | 976 | |
977 | if(mesonPdg == 111) SetTagBit(tag,kMCPi0); | |
36769d30 | 978 | else SetTagBit(tag,kMCEta); |
3c75ddf7 | 979 | } |
980 | ||
902aa95c | 981 | } |
982 | ||
8a2dbbff | 983 | //________________________________________________________________________________________________________ |
984 | void AliMCAnalysisUtils::CheckOverlapped2GammaDecay(const Int_t *labels, Int_t nlabels, Int_t mesonIndex, | |
985 | const TClonesArray *mcparticles, Int_t & tag ) | |
840124fe | 986 | { |
0d360687 | 987 | // Check if cluster is formed from the contribution of 2 decay photons from pi0 or eta. Input in AODMCParticles. |
f8006433 | 988 | |
2644ead9 | 989 | if(labels[0] < 0 || labels[0] > mcparticles->GetEntriesFast() || nlabels <= 1) |
990 | { | |
3c75ddf7 | 991 | if(fDebug > 2) |
992 | printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : label[0] %d, n primaries %d, nlabels %d \n", | |
993 | labels[0],mcparticles->GetEntriesFast(), nlabels); | |
994 | return; | |
995 | } | |
c5693f62 | 996 | |
f8006433 | 997 | AliAODMCParticle * meson = (AliAODMCParticle *) mcparticles->At(mesonIndex); |
3c75ddf7 | 998 | Int_t mesonPdg = meson->GetPdgCode(); |
2644ead9 | 999 | if(mesonPdg != 111 && mesonPdg != 221) |
1000 | { | |
3c75ddf7 | 1001 | printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Wrong pi0/eta PDG : %d \n",mesonPdg); |
1002 | return; | |
1003 | } | |
1004 | ||
1005 | if(fDebug > 2) printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - pdg %d, label %d, ndaughters %d\n", mesonPdg, mesonIndex, meson->GetNDaughters()); | |
1006 | ||
1007 | ||
1008 | //Get the daughters | |
2644ead9 | 1009 | if(meson->GetNDaughters() != 2) |
1010 | { | |
3c75ddf7 | 1011 | if(fDebug > 2) |
36769d30 | 1012 | printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(ESD) - Not overalapped. Number of daughters is %d, not 2 \n",meson->GetNDaughters()); |
3c75ddf7 | 1013 | return; |
1014 | } | |
2644ead9 | 1015 | |
3c75ddf7 | 1016 | Int_t iPhoton0 = meson->GetDaughter(0); |
1017 | Int_t iPhoton1 = meson->GetDaughter(1); | |
1018 | //if((iPhoton0 == -1) || (iPhoton1 == -1)){ | |
1019 | // if(fDebug > 2) | |
1020 | // printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Exit : Not overlapped. At least a daughter do not exists : d1 %d, d2 %d \n", iPhoton0, iPhoton1); | |
1021 | // return; | |
1022 | //} | |
1023 | AliAODMCParticle *photon0 = (AliAODMCParticle *) mcparticles->At(iPhoton0); | |
1024 | AliAODMCParticle *photon1 = (AliAODMCParticle *) mcparticles->At(iPhoton1); | |
2644ead9 | 1025 | |
3c75ddf7 | 1026 | //Check if both daughters are photons |
2644ead9 | 1027 | if(photon0->GetPdgCode() != 22 && photon1->GetPdgCode()!=22) |
1028 | { | |
3c75ddf7 | 1029 | printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Not overlapped. PDG: daughter 1 = %d, of daughter 2 = %d \n",photon0->GetPdgCode(),photon1->GetPdgCode()); |
1030 | return; | |
1031 | } | |
1032 | ||
1033 | if(fDebug > 2) | |
1034 | printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Daughter labels : photon0 = %d, photon1 = %d \n",iPhoton0,iPhoton1); | |
1035 | ||
1036 | //Check if both photons contribute to the cluster | |
1037 | Bool_t okPhoton0 = kFALSE; | |
1038 | Bool_t okPhoton1 = kFALSE; | |
1039 | ||
1040 | if(fDebug > 3) | |
1041 | printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Labels loop:\n"); | |
1042 | ||
36769d30 | 1043 | Bool_t conversion = kFALSE; |
1044 | ||
2644ead9 | 1045 | for(Int_t i = 0; i < nlabels; i++) |
1046 | { | |
1047 | if(fDebug > 3) | |
3c75ddf7 | 1048 | printf("\t label %d/%d: %d, ok? %d, %d\n", i, nlabels, labels[i], okPhoton0, okPhoton1); |
f8006433 | 1049 | |
2644ead9 | 1050 | if(labels[i]<0) continue; |
1051 | ||
3c75ddf7 | 1052 | //If we already found both, break the loop |
1053 | if(okPhoton0 && okPhoton1) break; | |
f8006433 | 1054 | |
3c75ddf7 | 1055 | Int_t index = labels[i]; |
36769d30 | 1056 | if (iPhoton0 == index) |
1057 | { | |
3c75ddf7 | 1058 | okPhoton0 = kTRUE; |
1059 | continue; | |
1060 | } | |
36769d30 | 1061 | else if (iPhoton1 == index) |
1062 | { | |
3c75ddf7 | 1063 | okPhoton1 = kTRUE; |
1064 | continue; | |
1065 | } | |
1066 | ||
1067 | //Trace back the mother in case it was a conversion | |
1068 | ||
36769d30 | 1069 | if(index >= mcparticles->GetEntriesFast()) |
1070 | { | |
3c75ddf7 | 1071 | printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) Particle index %d larger than size of list %d !!\n",index,mcparticles->GetEntriesFast()); |
1072 | continue; | |
1073 | } | |
1074 | ||
1075 | AliAODMCParticle * daught = (AliAODMCParticle*) mcparticles->At(index); | |
1076 | Int_t tmpindex = daught->GetMother(); | |
1077 | if(fDebug > 3) | |
1078 | printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Conversion? : mother %d\n",tmpindex); | |
1079 | ||
1080 | while(tmpindex>=0){ | |
1081 | ||
1082 | //MC particle of interest is the mother | |
1083 | if(fDebug > 3) | |
1084 | printf("\t parent index %d\n",tmpindex); | |
1085 | daught = (AliAODMCParticle*) mcparticles->At(tmpindex); | |
1086 | //printf("tmpindex %d\n",tmpindex); | |
36769d30 | 1087 | if (iPhoton0 == tmpindex) |
1088 | { | |
1089 | conversion = kTRUE; | |
1090 | okPhoton0 = kTRUE; | |
3c75ddf7 | 1091 | break; |
1092 | } | |
36769d30 | 1093 | else if (iPhoton1 == tmpindex) |
1094 | { | |
1095 | conversion = kTRUE; | |
1096 | okPhoton1 = kTRUE; | |
3c75ddf7 | 1097 | break; |
36769d30 | 1098 | } |
1099 | ||
3c75ddf7 | 1100 | tmpindex = daught->GetMother(); |
36769d30 | 1101 | |
3c75ddf7 | 1102 | }//While to check if pi0/eta daughter was one of these contributors to the cluster |
1103 | ||
2644ead9 | 1104 | if(i == 0 && (!okPhoton0 && !okPhoton1) && fDebug>=0 ) |
3c75ddf7 | 1105 | printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - Something happens, first label should be from a photon decay!\n"); |
1106 | ||
1107 | }//loop on list of labels | |
1108 | ||
1109 | //If both photons contribute tag as the corresponding meson. | |
36769d30 | 1110 | if(okPhoton0 && okPhoton1) |
1111 | { | |
3c75ddf7 | 1112 | if(fDebug > 2) |
36769d30 | 1113 | printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) - %s OVERLAPPED DECAY \n", |
1114 | (TDatabasePDG::Instance()->GetParticle(mesonPdg))->GetName()); | |
1115 | ||
1116 | if(!CheckTagBit(tag,kMCConversion) && conversion) | |
1117 | { | |
1118 | if(fDebug > 2) | |
1119 | printf("AliMCAnalysisUtils::CheckOverlapped2GammaDecay(AOD) Second decay photon produced a conversion \n"); | |
1120 | SetTagBit(tag,kMCConversion) ; | |
1121 | } | |
3c75ddf7 | 1122 | |
1123 | if(mesonPdg == 111) SetTagBit(tag,kMCPi0); | |
1124 | else SetTagBit(tag,kMCEta); | |
1125 | } | |
1126 | ||
902aa95c | 1127 | } |
7cd4e982 | 1128 | |
3e57baab | 1129 | //______________________________________________________________________________________________________ |
1130 | void AliMCAnalysisUtils::CheckLostDecayPair(const TObjArray* arrayCluster, Int_t iMom, Int_t iParent, | |
1131 | AliStack * stack, Int_t & tag) | |
1132 | { | |
0d360687 | 1133 | // Check on ESDs if the current decay photon has the second photon companion lost. |
3e57baab | 1134 | |
1135 | if(!arrayCluster || iMom < 0 || iParent < 0|| !stack) return; | |
1136 | ||
1137 | TParticle * parent= stack->Particle(iParent); | |
1138 | ||
1139 | if(parent->GetNDaughters()!=2) | |
1140 | { | |
1141 | SetTagBit(tag, kMCDecayPairLost); | |
1142 | return ; | |
1143 | } | |
1144 | ||
1145 | Int_t pairLabel = -1; | |
1146 | if ( iMom != parent->GetDaughter(0) ) pairLabel = parent->GetDaughter(0); | |
1147 | else if( iMom != parent->GetDaughter(1) ) pairLabel = parent->GetDaughter(1); | |
1148 | ||
1149 | if(pairLabel<0) | |
1150 | { | |
1151 | SetTagBit(tag, kMCDecayPairLost); | |
1152 | return ; | |
1153 | } | |
1154 | ||
1155 | for(Int_t iclus = 0; iclus < arrayCluster->GetEntriesFast(); iclus++) | |
1156 | { | |
1157 | AliVCluster * cluster = (AliVCluster*) arrayCluster->At(iclus); | |
1158 | for(UInt_t ilab = 0; ilab< cluster->GetNLabels(); ilab++) | |
1159 | { | |
1160 | Int_t label = cluster->GetLabels()[ilab]; | |
1161 | if(label==pairLabel) | |
1162 | { | |
1163 | SetTagBit(tag, kMCDecayPairInCalo); | |
1164 | return ; | |
1165 | } | |
1166 | else if(label== iParent || label== iMom) | |
1167 | { | |
1168 | continue; | |
1169 | } | |
1170 | else // check the ancestry | |
1171 | { | |
1172 | TParticle * mother = stack->Particle(label); | |
1173 | Int_t momPDG = TMath::Abs(mother->GetPdgCode()); | |
1174 | if(momPDG!=11 && momPDG!=22) continue; | |
1175 | ||
1176 | //Check if "mother" of entity is converted, if not, get the first non converted mother | |
1177 | Int_t iParentClus = mother->GetFirstMother(); | |
1178 | if(iParentClus < 0) continue; | |
1179 | ||
1180 | TParticle * parentClus = stack->Particle(iParentClus); | |
1181 | if(!parentClus) continue; | |
1182 | ||
1183 | Int_t parentClusPDG = TMath::Abs(parentClus->GetPdgCode()); | |
1184 | Int_t parentClusStatus = parentClus->GetStatusCode(); | |
1185 | ||
1186 | if( parentClusPDG != 22 && parentClusPDG != 11 && parentClusStatus != 0) continue; | |
1187 | ||
1188 | //printf("Conversion\n"); | |
1189 | ||
1190 | //Check if the mother is photon or electron with status not stable | |
1191 | while ((parentClusPDG == 22 || parentClusPDG == 11) && parentClusStatus != 1) | |
1192 | { | |
1193 | //New Mother | |
1194 | label = iParentClus; | |
1195 | momPDG = parentClusPDG; | |
1196 | ||
1197 | iParentClus = parentClus->GetFirstMother(); | |
1198 | if(iParentClus < 0) break; | |
1199 | ||
1200 | parentClus = stack->Particle(iParentClus); | |
1201 | if(!parentClus) break; | |
1202 | ||
1203 | parentClusPDG = TMath::Abs(parentClus->GetPdgCode()); | |
1204 | parentClusStatus = parentClus->GetStatusCode() ; | |
1205 | }//while | |
1206 | ||
1207 | if((momPDG == 22 || parentClusPDG ==22) && (label==pairLabel || iParentClus == pairLabel)) | |
1208 | { | |
1209 | SetTagBit(tag, kMCDecayPairInCalo); | |
1210 | //printf("Conversion is paired\n"); | |
1211 | return ; | |
1212 | } | |
1213 | else continue; | |
1214 | ||
1215 | } | |
1216 | } | |
1217 | } // cluster loop | |
1218 | ||
1219 | SetTagBit(tag, kMCDecayPairLost); | |
1220 | ||
1221 | } | |
1222 | ||
1223 | //______________________________________________________________________________________________________ | |
1224 | void AliMCAnalysisUtils::CheckLostDecayPair(const TObjArray * arrayCluster,Int_t iMom, Int_t iParent, | |
1225 | const TClonesArray* mcparticles, Int_t & tag) | |
1226 | { | |
0d360687 | 1227 | // Check on AODs if the current decay photon has the second photon companion lost. |
3e57baab | 1228 | |
1229 | if(!arrayCluster || iMom < 0 || iParent < 0|| !mcparticles) return; | |
1230 | ||
1231 | AliAODMCParticle * parent = (AliAODMCParticle*) mcparticles->At(iParent); | |
1232 | ||
1233 | //printf("*** Check label %d with parent %d\n",iMom, iParent); | |
1234 | ||
1235 | if(parent->GetNDaughters()!=2) | |
1236 | { | |
1237 | SetTagBit(tag, kMCDecayPairLost); | |
1238 | //printf("\t ndaugh = %d\n",parent->GetNDaughters()); | |
1239 | return ; | |
1240 | } | |
1241 | ||
1242 | Int_t pairLabel = -1; | |
1243 | if ( iMom != parent->GetDaughter(0) ) pairLabel = parent->GetDaughter(0); | |
1244 | else if( iMom != parent->GetDaughter(1) ) pairLabel = parent->GetDaughter(1); | |
1245 | ||
1246 | if(pairLabel<0) | |
1247 | { | |
1248 | //printf("\t pair Label not found = %d\n",pairLabel); | |
1249 | SetTagBit(tag, kMCDecayPairLost); | |
1250 | return ; | |
1251 | } | |
1252 | ||
1253 | //printf("\t *** find pair %d\n",pairLabel); | |
1254 | ||
1255 | for(Int_t iclus = 0; iclus < arrayCluster->GetEntriesFast(); iclus++) | |
1256 | { | |
1257 | AliVCluster * cluster = (AliVCluster*) arrayCluster->At(iclus); | |
1258 | //printf("\t \t ** Cluster %d, nlabels %d\n",iclus,cluster->GetNLabels()); | |
1259 | for(UInt_t ilab = 0; ilab< cluster->GetNLabels(); ilab++) | |
1260 | { | |
1261 | Int_t label = cluster->GetLabels()[ilab]; | |
1262 | ||
1263 | //printf("\t \t label %d\n",label); | |
1264 | ||
1265 | if(label==pairLabel) | |
1266 | { | |
1267 | //printf("\t \t Pair found\n"); | |
1268 | SetTagBit(tag, kMCDecayPairInCalo); | |
1269 | return ; | |
1270 | } | |
1271 | else if(label== iParent || label== iMom) | |
1272 | { | |
1273 | //printf("\t \t skip\n"); | |
1274 | continue; | |
1275 | } | |
1276 | else // check the ancestry | |
1277 | { | |
1278 | AliAODMCParticle * mother = (AliAODMCParticle*) mcparticles->At(label); | |
1279 | Int_t momPDG = TMath::Abs(mother->GetPdgCode()); | |
1280 | if(momPDG!=11 && momPDG!=22) | |
1281 | { | |
1282 | //printf("\t \t skip, pdg %d\n",momPDG); | |
1283 | continue; | |
1284 | } | |
1285 | ||
1286 | //Check if "mother" of entity is converted, if not, get the first non converted mother | |
1287 | Int_t iParentClus = mother->GetMother(); | |
1288 | if(iParentClus < 0) continue; | |
1289 | ||
1290 | AliAODMCParticle * parentClus = (AliAODMCParticle*) mcparticles->At(iParentClus); | |
1291 | if(!parentClus) continue; | |
1292 | ||
1293 | Int_t parentClusPDG = TMath::Abs(parentClus->GetPdgCode()); | |
1294 | Int_t parentClusStatus = parentClus->GetStatus(); | |
1295 | ||
1296 | if( parentClusPDG != 22 && parentClusPDG != 11 && parentClusStatus != 0) | |
1297 | { | |
1298 | //printf("\t \t skip, not a conversion, parent: pdg %d, status %d\n",parentClusPDG,parentClusStatus); | |
1299 | continue; | |
1300 | } | |
1301 | ||
1302 | //printf("\t \t Conversion\n"); | |
1303 | ||
1304 | //Check if the mother is photon or electron with status not stable | |
1305 | while ((parentClusPDG == 22 || parentClusPDG == 11) && parentClusStatus != 1) | |
1306 | { | |
1307 | //New Mother | |
1308 | label = iParentClus; | |
1309 | momPDG = parentClusPDG; | |
1310 | ||
1311 | iParentClus = parentClus->GetMother(); | |
1312 | if(iParentClus < 0) break; | |
1313 | ||
1314 | parentClus = (AliAODMCParticle*) mcparticles->At(iParentClus); | |
1315 | if(!parentClus) break; | |
1316 | ||
1317 | parentClusPDG = TMath::Abs(parentClus->GetPdgCode()); | |
1318 | parentClusStatus = parentClus->GetStatus() ; | |
1319 | }//while | |
1320 | ||
1321 | if((momPDG == 22 || parentClusPDG ==22) && (label==pairLabel || iParentClus == pairLabel)) | |
1322 | { | |
1323 | SetTagBit(tag, kMCDecayPairInCalo); | |
1324 | //printf("\t \t Conversion is paired: mom %d, parent %d\n",label,iParentClus); | |
1325 | return ; | |
1326 | } | |
1327 | else | |
1328 | { | |
1329 | //printf("\t \t Skip, finally label %d, pdg %d, parent label %d, pdg %d, status %d\n",label,momPDG,iParentClus,parentClusPDG,parentClusStatus); | |
1330 | continue; | |
1331 | } | |
1332 | ||
1333 | } | |
1334 | } | |
1335 | } // cluster loop | |
1336 | ||
1337 | ||
1338 | SetTagBit(tag, kMCDecayPairLost); | |
1339 | ||
1340 | ||
1341 | } | |
1342 | ||
8a2dbbff | 1343 | //_____________________________________________________________________ |
c5693f62 | 1344 | TList * AliMCAnalysisUtils::GetJets(const AliCaloTrackReader * reader) |
840124fe | 1345 | { |
0d360687 | 1346 | // Return list of jets (TParticles) and index of most likely parton that originated it. |
1347 | ||
7cd4e982 | 1348 | AliStack * stack = reader->GetStack(); |
1349 | Int_t iEvent = reader->GetEventNumber(); | |
1350 | AliGenEventHeader * geh = reader->GetGenEventHeader(); | |
1351 | if(fCurrentEvent!=iEvent){ | |
1352 | fCurrentEvent = iEvent; | |
1353 | fJetsList = new TList; | |
1354 | Int_t nTriggerJets = 0; | |
1355 | Float_t tmpjet[]={0,0,0,0}; | |
1356 | ||
1357 | //printf("Event %d %d\n",fCurrentEvent,iEvent); | |
1358 | //Get outgoing partons | |
1359 | if(stack->GetNtrack() < 8) return fJetsList; | |
1360 | TParticle * parton1 = stack->Particle(6); | |
1361 | TParticle * parton2 = stack->Particle(7); | |
1362 | if(fDebug > 2){ | |
1363 | printf("AliMCAnalysisUtils::GetJets() - parton 6 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n", | |
f8006433 | 1364 | parton1->GetName(),parton1->Pt(),parton1->Energy(),parton1->Phi()*TMath::RadToDeg(),parton1->Eta()); |
7cd4e982 | 1365 | printf("AliMCAnalysisUtils::GetJets() - parton 7 : %s, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n", |
f8006433 | 1366 | parton2->GetName(),parton2->Pt(),parton2->Energy(),parton2->Phi()*TMath::RadToDeg(),parton2->Eta()); |
7cd4e982 | 1367 | } |
f8006433 | 1368 | // //Trace the jet from the mother parton |
1369 | // Float_t pt = 0; | |
1370 | // Float_t pt1 = 0; | |
1371 | // Float_t pt2 = 0; | |
1372 | // Float_t e = 0; | |
1373 | // Float_t e1 = 0; | |
1374 | // Float_t e2 = 0; | |
1375 | // TParticle * tmptmp = new TParticle; | |
1376 | // for(Int_t i = 0; i< stack->GetNprimary(); i++){ | |
1377 | // tmptmp = stack->Particle(i); | |
7cd4e982 | 1378 | |
f8006433 | 1379 | // if(tmptmp->GetStatusCode() == 1){ |
1380 | // pt = tmptmp->Pt(); | |
1381 | // e = tmptmp->Energy(); | |
1382 | // Int_t imom = tmptmp->GetFirstMother(); | |
1383 | // Int_t imom1 = 0; | |
1384 | // //printf("1st imom %d\n",imom); | |
1385 | // while(imom > 5){ | |
1386 | // imom1=imom; | |
1387 | // tmptmp = stack->Particle(imom); | |
1388 | // imom = tmptmp->GetFirstMother(); | |
1389 | // //printf("imom %d \n",imom); | |
1390 | // } | |
1391 | // //printf("Last imom %d %d\n",imom1, imom); | |
1392 | // if(imom1 == 6) { | |
1393 | // pt1+=pt; | |
1394 | // e1+=e; | |
1395 | // } | |
1396 | // else if (imom1 == 7){ | |
1397 | // pt2+=pt; | |
1398 | // e2+=e; } | |
1399 | // }// status 1 | |
1400 | ||
1401 | // }// for | |
7cd4e982 | 1402 | |
f8006433 | 1403 | // printf("JET 1, pt %2.2f, e %2.2f; JET 2, pt %2.2f, e %2.2f \n",pt1,e1,pt2,e2); |
7cd4e982 | 1404 | |
1405 | //Get the jet, different way for different generator | |
1406 | //PYTHIA | |
0d360687 | 1407 | if(fMCGenerator == kPythia) |
4d52b003 | 1408 | { |
f8006433 | 1409 | TParticle * jet = 0x0; |
7cd4e982 | 1410 | AliGenPythiaEventHeader* pygeh= (AliGenPythiaEventHeader*) geh; |
1411 | nTriggerJets = pygeh->NTriggerJets(); | |
1412 | if(fDebug > 1) | |
f8006433 | 1413 | printf("AliMCAnalysisUtils::GetJets() - PythiaEventHeader: Njets: %d\n",nTriggerJets); |
1414 | ||
265cab75 | 1415 | for(Int_t i = 0; i< nTriggerJets; i++) |
1416 | { | |
f8006433 | 1417 | pygeh->TriggerJet(i, tmpjet); |
1418 | jet = new TParticle(94, 21, -1, -1, -1, -1, tmpjet[0],tmpjet[1],tmpjet[2],tmpjet[3], 0,0,0,0); | |
1419 | //Assign an outgoing parton as mother | |
1420 | Float_t phidiff1 = TMath::Abs(jet->Phi()-parton1->Phi()); | |
1421 | Float_t phidiff2 = TMath::Abs(jet->Phi()-parton2->Phi()); | |
1422 | if(phidiff1 > phidiff2) jet->SetFirstMother(7); | |
1423 | else jet->SetFirstMother(6); | |
1424 | //jet->Print(); | |
1425 | if(fDebug > 1) | |
1426 | printf("AliMCAnalysisUtils::GetJets() - PYTHIA Jet %d: mother %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n", | |
1427 | i, jet->GetFirstMother(),jet->Pt(),jet->Energy(),jet->Phi()*TMath::RadToDeg(),jet->Eta()); | |
1428 | fJetsList->Add(jet); | |
7cd4e982 | 1429 | } |
1430 | }//Pythia triggered jets | |
1431 | //HERWIG | |
0d360687 | 1432 | else if (fMCGenerator == kHerwig) |
265cab75 | 1433 | { |
7cd4e982 | 1434 | Int_t pdg = -1; |
1435 | //Check parton 1 | |
1436 | TParticle * tmp = parton1; | |
265cab75 | 1437 | if(parton1->GetPdgCode()!=22) |
1438 | { | |
f8006433 | 1439 | while(pdg != 94){ |
1440 | if(tmp->GetFirstDaughter()==-1) return fJetsList; | |
1441 | tmp = stack->Particle(tmp->GetFirstDaughter()); | |
1442 | pdg = tmp->GetPdgCode(); | |
1443 | }//while | |
1444 | ||
1445 | //Add found jet to list | |
1446 | TParticle *jet1 = new TParticle(*tmp); | |
1447 | jet1->SetFirstMother(6); | |
1448 | fJetsList->Add(jet1); | |
1449 | //printf("jet 1: first daughter %d, last daughter %d\n", tmp->GetFirstDaughter(), tmp->GetLastDaughter()); | |
1450 | //tmp = stack->Particle(tmp->GetFirstDaughter()); | |
1451 | //tmp->Print(); | |
1452 | //jet1->Print(); | |
1453 | if(fDebug > 1) | |
1454 | printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 1: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n", | |
1455 | jet1->GetFirstMother(),jet1->GetStatusCode(),jet1->Pt(),jet1->Energy(),jet1->Phi()*TMath::RadToDeg(),jet1->Eta()); | |
7cd4e982 | 1456 | }//not photon |
1457 | ||
1458 | //Check parton 2 | |
1459 | pdg = -1; | |
1460 | tmp = parton2; | |
265cab75 | 1461 | if(parton2->GetPdgCode()!=22) |
1462 | { | |
1463 | while(pdg != 94) | |
1464 | { | |
f8006433 | 1465 | if(tmp->GetFirstDaughter()==-1) return fJetsList; |
f8006433 | 1466 | tmp = stack->Particle(tmp->GetFirstDaughter()); |
1467 | pdg = tmp->GetPdgCode(); | |
1468 | }//while | |
265cab75 | 1469 | |
f8006433 | 1470 | //Add found jet to list |
1471 | TParticle *jet2 = new TParticle(*tmp); | |
1472 | jet2->SetFirstMother(7); | |
1473 | fJetsList->Add(jet2); | |
1474 | //jet2->Print(); | |
1475 | if(fDebug > 1) | |
1476 | printf("AliMCAnalysisUtils::GetJets() - HERWIG Jet 2: mother %d, status %d, pt %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n", | |
1477 | jet2->GetFirstMother(),jet2->GetStatusCode(),jet2->Pt(),jet2->Energy(),jet2->Phi()*TMath::RadToDeg(),jet2->Eta()); | |
1478 | //Int_t first = tmp->GetFirstDaughter(); | |
1479 | //Int_t last = tmp->GetLastDaughter(); | |
1480 | //printf("jet 2: first daughter %d, last daughter %d, pdg %d\n",first, last, tmp->GetPdgCode()); | |
7cd4e982 | 1481 | // for(Int_t d = first ; d < last+1; d++){ |
f8006433 | 1482 | // tmp = stack->Particle(d); |
1483 | // if(i == tmp->GetFirstMother()) | |
1484 | // printf("Daughter n %d, Mother %d, name %s, status %d, pT %2.2f,E %2.2f, phi %2.2f, eta %2.2f \n", | |
1485 | // d,tmp->GetFirstMother(), tmp->GetName(), tmp->GetStatusCode(),tmp->Pt(),tmp->Energy(),tmp->Phi()*TMath::RadToDeg(),tmp->Eta()); | |
1486 | // } | |
1487 | //tmp->Print(); | |
7cd4e982 | 1488 | }//not photon |
1489 | }//Herwig generated jets | |
1490 | } | |
1491 | ||
1492 | return fJetsList; | |
1493 | } | |
1494 | ||
51a0ace5 | 1495 | |
8a2dbbff | 1496 | //__________________________________________________________________________________________________________ |
1497 | TLorentzVector AliMCAnalysisUtils::GetDaughter(Int_t idaugh, Int_t label, | |
64be16b4 | 1498 | const AliCaloTrackReader* reader, |
1499 | Int_t & pdg, Int_t & status, Bool_t & ok, Int_t & daughlabel) | |
1500 | { | |
0d360687 | 1501 | // Return the kinematics of the particle that generated the signal, its pdg and its status and its label mother. |
1502 | fDaughMom.SetPxPyPzE(0,0,0,0); | |
64be16b4 | 1503 | |
1504 | if(reader->ReadStack()) | |
1505 | { | |
1506 | if(!reader->GetStack()) | |
1507 | { | |
1508 | if (fDebug >=0) | |
1509 | printf("AliMCAnalysisUtils::GetDaughter() - Stack is not available, check analysis settings in configuration file, STOP!!\n"); | |
1510 | ||
1511 | ok=kFALSE; | |
0d360687 | 1512 | return fDaughMom; |
64be16b4 | 1513 | } |
4d6b53b6 | 1514 | |
1515 | Int_t nprimaries = reader->GetStack()->GetNtrack(); | |
1516 | if(label >= 0 && label < nprimaries) | |
64be16b4 | 1517 | { |
1518 | TParticle * momP = reader->GetStack()->Particle(label); | |
1519 | daughlabel = momP->GetDaughter(idaugh); | |
1520 | ||
4d6b53b6 | 1521 | if(daughlabel < 0 || daughlabel >= nprimaries) |
1522 | { | |
1523 | ok = kFALSE; | |
0d360687 | 1524 | return fDaughMom; |
4d6b53b6 | 1525 | } |
1526 | ||
64be16b4 | 1527 | TParticle * daughP = reader->GetStack()->Particle(daughlabel); |
0d360687 | 1528 | daughP->Momentum(fDaughMom); |
64be16b4 | 1529 | pdg = daughP->GetPdgCode(); |
1530 | status = daughP->GetStatusCode(); | |
1531 | } | |
1532 | else | |
1533 | { | |
1534 | ok = kFALSE; | |
0d360687 | 1535 | return fDaughMom; |
64be16b4 | 1536 | } |
1537 | } | |
1538 | else if(reader->ReadAODMCParticles()) | |
1539 | { | |
1540 | TClonesArray* mcparticles = reader->GetAODMCParticles(); | |
1541 | if(!mcparticles) | |
1542 | { | |
1543 | if(fDebug >= 0) | |
1544 | printf("AliMCAnalysisUtils::GetDaughter() - AODMCParticles is not available, check analysis settings in configuration file!!\n"); | |
1545 | ||
1546 | ok=kFALSE; | |
0d360687 | 1547 | return fDaughMom; |
64be16b4 | 1548 | } |
1549 | ||
1550 | Int_t nprimaries = mcparticles->GetEntriesFast(); | |
1551 | if(label >= 0 && label < nprimaries) | |
1552 | { | |
1553 | AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label); | |
1554 | daughlabel = momP->GetDaughter(idaugh); | |
1555 | ||
4d6b53b6 | 1556 | if(daughlabel < 0 || daughlabel >= nprimaries) |
1557 | { | |
1558 | ok = kFALSE; | |
0d360687 | 1559 | return fDaughMom; |
4d6b53b6 | 1560 | } |
1561 | ||
64be16b4 | 1562 | AliAODMCParticle * daughP = (AliAODMCParticle *) mcparticles->At(daughlabel); |
0d360687 | 1563 | fDaughMom.SetPxPyPzE(daughP->Px(),daughP->Py(),daughP->Pz(),daughP->E()); |
64be16b4 | 1564 | pdg = daughP->GetPdgCode(); |
1565 | status = daughP->GetStatus(); | |
1566 | } | |
1567 | else | |
1568 | { | |
1569 | ok = kFALSE; | |
0d360687 | 1570 | return fDaughMom; |
64be16b4 | 1571 | } |
1572 | } | |
1573 | ||
1574 | ok = kTRUE; | |
1575 | ||
0d360687 | 1576 | return fDaughMom; |
64be16b4 | 1577 | } |
1578 | ||
8a2dbbff | 1579 | //______________________________________________________________________________________________________ |
1580 | TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader, Bool_t & ok) | |
9036b2a6 | 1581 | { |
0d360687 | 1582 | // Return the kinematics of the particle that generated the signal. |
9036b2a6 | 1583 | |
64be16b4 | 1584 | Int_t pdg = -1; Int_t status = -1; Int_t momlabel = -1; |
1585 | return GetMother(label,reader,pdg,status, ok,momlabel); | |
1586 | } | |
1587 | ||
8a2dbbff | 1588 | //_________________________________________________________________________________________ |
1589 | TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader, | |
64be16b4 | 1590 | Int_t & pdg, Int_t & status, Bool_t & ok) |
1591 | { | |
0d360687 | 1592 | // Return the kinematics of the particle that generated the signal. |
64be16b4 | 1593 | |
1594 | Int_t momlabel = -1; | |
1595 | return GetMother(label,reader,pdg,status, ok,momlabel); | |
51a0ace5 | 1596 | } |
1597 | ||
8a2dbbff | 1598 | //______________________________________________________________________________________________________ |
1599 | TLorentzVector AliMCAnalysisUtils::GetMother(Int_t label, const AliCaloTrackReader* reader, | |
64be16b4 | 1600 | Int_t & pdg, Int_t & status, Bool_t & ok, Int_t & momlabel) |
51a0ace5 | 1601 | { |
0d360687 | 1602 | // Return the kinematics of the particle that generated the signal, its pdg and its status and its label mother. |
51a0ace5 | 1603 | |
0d360687 | 1604 | fMotherMom.SetPxPyPzE(0,0,0,0); |
9036b2a6 | 1605 | |
1606 | if(reader->ReadStack()) | |
1607 | { | |
51a0ace5 | 1608 | if(!reader->GetStack()) |
1609 | { | |
9036b2a6 | 1610 | if (fDebug >=0) |
1611 | printf("AliMCAnalysisUtils::GetMother() - Stack is not available, check analysis settings in configuration file, STOP!!\n"); | |
51a0ace5 | 1612 | |
1613 | ok=kFALSE; | |
0d360687 | 1614 | return fMotherMom; |
9036b2a6 | 1615 | } |
1616 | if(label >= 0 && label < reader->GetStack()->GetNtrack()) | |
1617 | { | |
1618 | TParticle * momP = reader->GetStack()->Particle(label); | |
0d360687 | 1619 | momP->Momentum(fMotherMom); |
1620 | pdg = momP->GetPdgCode(); | |
1621 | status = momP->GetStatusCode(); | |
64be16b4 | 1622 | momlabel = momP->GetFirstMother(); |
51a0ace5 | 1623 | } |
1624 | else | |
1625 | { | |
1626 | ok = kFALSE; | |
0d360687 | 1627 | return fMotherMom; |
9036b2a6 | 1628 | } |
1629 | } | |
1630 | else if(reader->ReadAODMCParticles()) | |
1631 | { | |
2644ead9 | 1632 | TClonesArray* mcparticles = reader->GetAODMCParticles(); |
9036b2a6 | 1633 | if(!mcparticles) |
1634 | { | |
1635 | if(fDebug >= 0) | |
1636 | printf("AliMCAnalysisUtils::GetMother() - AODMCParticles is not available, check analysis settings in configuration file!!\n"); | |
51a0ace5 | 1637 | |
1638 | ok=kFALSE; | |
0d360687 | 1639 | return fMotherMom; |
9036b2a6 | 1640 | } |
1641 | ||
1642 | Int_t nprimaries = mcparticles->GetEntriesFast(); | |
1643 | if(label >= 0 && label < nprimaries) | |
1644 | { | |
1645 | AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label); | |
0d360687 | 1646 | fMotherMom.SetPxPyPzE(momP->Px(),momP->Py(),momP->Pz(),momP->E()); |
1647 | pdg = momP->GetPdgCode(); | |
1648 | status = momP->GetStatus(); | |
64be16b4 | 1649 | momlabel = momP->GetMother(); |
1650 | } | |
51a0ace5 | 1651 | else |
1652 | { | |
1653 | ok = kFALSE; | |
0d360687 | 1654 | return fMotherMom; |
9036b2a6 | 1655 | } |
1656 | } | |
1657 | ||
51a0ace5 | 1658 | ok = kTRUE; |
1659 | ||
0d360687 | 1660 | return fMotherMom; |
9036b2a6 | 1661 | } |
1662 | ||
8a2dbbff | 1663 | //___________________________________________________________________________________ |
1664 | TLorentzVector AliMCAnalysisUtils::GetMotherWithPDG(Int_t label, Int_t pdg, | |
1665 | const AliCaloTrackReader* reader, | |
36769d30 | 1666 | Bool_t & ok, Int_t & momlabel) |
9036b2a6 | 1667 | { |
0d360687 | 1668 | // Return the kinematics of the particle that generated the signal. |
9036b2a6 | 1669 | |
0d360687 | 1670 | fGMotherMom.SetPxPyPzE(0,0,0,0); |
9036b2a6 | 1671 | |
1672 | if(reader->ReadStack()) | |
1673 | { | |
1674 | if(!reader->GetStack()) | |
1675 | { | |
1676 | if (fDebug >=0) | |
1677 | printf("AliMCAnalysisUtils::GetMotherWithPDG() - Stack is not available, check analysis settings in configuration file, STOP!!\n"); | |
51a0ace5 | 1678 | |
1679 | ok = kFALSE; | |
0d360687 | 1680 | return fGMotherMom; |
9036b2a6 | 1681 | } |
36769d30 | 1682 | |
9036b2a6 | 1683 | if(label >= 0 && label < reader->GetStack()->GetNtrack()) |
1684 | { | |
1685 | TParticle * momP = reader->GetStack()->Particle(label); | |
1686 | ||
1687 | Int_t grandmomLabel = momP->GetFirstMother(); | |
1688 | Int_t grandmomPDG = -1; | |
1689 | TParticle * grandmomP = 0x0; | |
1690 | while (grandmomLabel >=0 ) | |
1691 | { | |
1692 | grandmomP = reader->GetStack()->Particle(grandmomLabel); | |
1693 | grandmomPDG = grandmomP->GetPdgCode(); | |
1694 | if(grandmomPDG==pdg) | |
1695 | { | |
36769d30 | 1696 | momlabel = grandmomLabel; |
0d360687 | 1697 | fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy()); |
9036b2a6 | 1698 | break; |
1699 | } | |
1700 | ||
1701 | grandmomLabel = grandmomP->GetFirstMother(); | |
1702 | ||
1703 | } | |
1704 | ||
1705 | if(grandmomPDG!=pdg) printf("AliMCAnalysisUtils::GetMotherWithPDG(ESD) - mother with PDG %d, not found! \n",pdg); | |
1706 | } | |
1707 | } | |
1708 | else if(reader->ReadAODMCParticles()) | |
1709 | { | |
2644ead9 | 1710 | TClonesArray* mcparticles = reader->GetAODMCParticles(); |
9036b2a6 | 1711 | if(!mcparticles) |
1712 | { | |
1713 | if(fDebug >= 0) | |
1714 | printf("AliMCAnalysisUtils::GetMotherWithPDG() - AODMCParticles is not available, check analysis settings in configuration file!!\n"); | |
51a0ace5 | 1715 | |
1716 | ok=kFALSE; | |
0d360687 | 1717 | return fGMotherMom; |
9036b2a6 | 1718 | } |
1719 | ||
1720 | Int_t nprimaries = mcparticles->GetEntriesFast(); | |
1721 | if(label >= 0 && label < nprimaries) | |
1722 | { | |
1723 | AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label); | |
1724 | ||
1725 | Int_t grandmomLabel = momP->GetMother(); | |
1726 | Int_t grandmomPDG = -1; | |
1727 | AliAODMCParticle * grandmomP = 0x0; | |
1728 | while (grandmomLabel >=0 ) | |
1729 | { | |
1730 | grandmomP = (AliAODMCParticle *) mcparticles->At(grandmomLabel); | |
1731 | grandmomPDG = grandmomP->GetPdgCode(); | |
1732 | if(grandmomPDG==pdg) | |
1733 | { | |
1734 | //printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d FOUND! \n",pdg); | |
36769d30 | 1735 | momlabel = grandmomLabel; |
0d360687 | 1736 | fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E()); |
9036b2a6 | 1737 | break; |
1738 | } | |
1739 | ||
1740 | grandmomLabel = grandmomP->GetMother(); | |
1741 | ||
1742 | } | |
1743 | ||
1744 | if(grandmomPDG!=pdg) printf("AliMCAnalysisUtils::GetMotherWithPDG(AOD) - mother with PDG %d, NOT found! \n",pdg); | |
1745 | ||
1746 | } | |
1747 | } | |
1748 | ||
51a0ace5 | 1749 | ok = kTRUE; |
9036b2a6 | 1750 | |
0d360687 | 1751 | return fGMotherMom; |
9036b2a6 | 1752 | } |
7cd4e982 | 1753 | |
8a2dbbff | 1754 | //______________________________________________________________________________________________ |
1755 | TLorentzVector AliMCAnalysisUtils::GetGrandMother(Int_t label, const AliCaloTrackReader* reader, | |
914b9fe7 | 1756 | Int_t & pdg, Int_t & status, Bool_t & ok, |
1757 | Int_t & grandMomLabel, Int_t & greatMomLabel) | |
1758 | { | |
0d360687 | 1759 | // Return the kinematics of the particle that generated the signal. |
914b9fe7 | 1760 | |
0d360687 | 1761 | fGMotherMom.SetPxPyPzE(0,0,0,0); |
914b9fe7 | 1762 | |
1763 | if(reader->ReadStack()) | |
1764 | { | |
1765 | if(!reader->GetStack()) | |
1766 | { | |
1767 | if (fDebug >=0) | |
64be16b4 | 1768 | printf("AliMCAnalysisUtils::GetGrandMother() - Stack is not available, check analysis settings in configuration file, STOP!!\n"); |
914b9fe7 | 1769 | |
1770 | ok = kFALSE; | |
0d360687 | 1771 | return fGMotherMom; |
914b9fe7 | 1772 | } |
1773 | ||
1774 | if(label >= 0 && label < reader->GetStack()->GetNtrack()) | |
1775 | { | |
1776 | TParticle * momP = reader->GetStack()->Particle(label); | |
1777 | ||
1778 | grandMomLabel = momP->GetFirstMother(); | |
1779 | ||
1780 | TParticle * grandmomP = 0x0; | |
1781 | ||
1782 | if (grandMomLabel >=0 ) | |
1783 | { | |
1784 | grandmomP = reader->GetStack()->Particle(grandMomLabel); | |
1785 | pdg = grandmomP->GetPdgCode(); | |
1786 | status = grandmomP->GetStatusCode(); | |
1787 | ||
0d360687 | 1788 | fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->Energy()); |
914b9fe7 | 1789 | greatMomLabel = grandmomP->GetFirstMother(); |
1790 | ||
1791 | } | |
1792 | } | |
1793 | } | |
1794 | else if(reader->ReadAODMCParticles()) | |
1795 | { | |
1796 | TClonesArray* mcparticles = reader->GetAODMCParticles(); | |
1797 | if(!mcparticles) | |
1798 | { | |
1799 | if(fDebug >= 0) | |
64be16b4 | 1800 | printf("AliMCAnalysisUtils::GetGrandMother() - AODMCParticles is not available, check analysis settings in configuration file!!\n"); |
914b9fe7 | 1801 | |
1802 | ok=kFALSE; | |
0d360687 | 1803 | return fGMotherMom; |
914b9fe7 | 1804 | } |
1805 | ||
1806 | Int_t nprimaries = mcparticles->GetEntriesFast(); | |
1807 | if(label >= 0 && label < nprimaries) | |
1808 | { | |
1809 | AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label); | |
1810 | ||
1811 | grandMomLabel = momP->GetMother(); | |
1812 | ||
1813 | AliAODMCParticle * grandmomP = 0x0; | |
1814 | ||
1815 | if(grandMomLabel >=0 ) | |
1816 | { | |
1817 | grandmomP = (AliAODMCParticle *) mcparticles->At(grandMomLabel); | |
1818 | pdg = grandmomP->GetPdgCode(); | |
1819 | status = grandmomP->GetStatus(); | |
1820 | ||
0d360687 | 1821 | fGMotherMom.SetPxPyPzE(grandmomP->Px(),grandmomP->Py(),grandmomP->Pz(),grandmomP->E()); |
914b9fe7 | 1822 | greatMomLabel = grandmomP->GetMother(); |
1823 | ||
1824 | } | |
1825 | } | |
1826 | } | |
1827 | ||
1828 | ok = kTRUE; | |
1829 | ||
0d360687 | 1830 | return fGMotherMom; |
914b9fe7 | 1831 | } |
1832 | ||
8a2dbbff | 1833 | //_______________________________________________________________________________________________________________ |
1834 | void AliMCAnalysisUtils::GetMCDecayAsymmetryAngleForPDG(Int_t label, Int_t pdg, const AliCaloTrackReader* reader, | |
4efb6957 | 1835 | Float_t & asym, Float_t & angle, Bool_t & ok) |
8e81c2cf | 1836 | { |
0d360687 | 1837 | // In case of an eta or pi0 decay into 2 photons, get the asymmetry in the energy of the photons. |
8e81c2cf | 1838 | |
1839 | if(reader->ReadStack()) | |
1840 | { | |
1841 | if(!reader->GetStack()) | |
1842 | { | |
1843 | if (fDebug >=0) | |
64be16b4 | 1844 | printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG() - Stack is not available, check analysis settings in configuration file, STOP!!\n"); |
8e81c2cf | 1845 | |
1846 | ok = kFALSE; | |
8e81c2cf | 1847 | } |
1848 | if(label >= 0 && label < reader->GetStack()->GetNtrack()) | |
1849 | { | |
1850 | TParticle * momP = reader->GetStack()->Particle(label); | |
1851 | ||
1852 | Int_t grandmomLabel = momP->GetFirstMother(); | |
1853 | Int_t grandmomPDG = -1; | |
1854 | TParticle * grandmomP = 0x0; | |
1855 | while (grandmomLabel >=0 ) | |
1856 | { | |
1857 | grandmomP = reader->GetStack()->Particle(grandmomLabel); | |
1858 | grandmomPDG = grandmomP->GetPdgCode(); | |
1859 | ||
1860 | if(grandmomPDG==pdg) break; | |
1861 | ||
1862 | grandmomLabel = grandmomP->GetFirstMother(); | |
1863 | ||
1864 | } | |
1865 | ||
1866 | if(grandmomPDG==pdg && grandmomP->GetNDaughters()==2) | |
1867 | { | |
1868 | TParticle * d1 = reader->GetStack()->Particle(grandmomP->GetDaughter(0)); | |
1869 | TParticle * d2 = reader->GetStack()->Particle(grandmomP->GetDaughter(1)); | |
1870 | if(d1->GetPdgCode() == 22 && d1->GetPdgCode() == 22) | |
1871 | { | |
1872 | asym = (d1->Energy()-d2->Energy())/grandmomP->Energy(); | |
0d360687 | 1873 | d1->Momentum(fDaughMom ); |
1874 | d2->Momentum(fDaughMom2); | |
1875 | angle = fDaughMom.Angle(fDaughMom2.Vect()); | |
8e81c2cf | 1876 | } |
1877 | } | |
1878 | else | |
1879 | { | |
1880 | ok=kFALSE; | |
64be16b4 | 1881 | printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG(ESD) - mother with PDG %d, not found! \n",pdg); |
8e81c2cf | 1882 | } |
1883 | ||
1884 | } // good label | |
1885 | } | |
1886 | else if(reader->ReadAODMCParticles()) | |
1887 | { | |
2644ead9 | 1888 | TClonesArray* mcparticles = reader->GetAODMCParticles(); |
8e81c2cf | 1889 | if(!mcparticles) |
1890 | { | |
1891 | if(fDebug >= 0) | |
64be16b4 | 1892 | printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG() - AODMCParticles is not available, check analysis settings in configuration file!!\n"); |
8e81c2cf | 1893 | |
1894 | ok=kFALSE; | |
c31af35c | 1895 | return; |
8e81c2cf | 1896 | } |
1897 | ||
1898 | Int_t nprimaries = mcparticles->GetEntriesFast(); | |
1899 | if(label >= 0 && label < nprimaries) | |
1900 | { | |
1901 | AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label); | |
1902 | ||
1903 | Int_t grandmomLabel = momP->GetMother(); | |
1904 | Int_t grandmomPDG = -1; | |
1905 | AliAODMCParticle * grandmomP = 0x0; | |
1906 | while (grandmomLabel >=0 ) | |
1907 | { | |
1908 | grandmomP = (AliAODMCParticle *) mcparticles->At(grandmomLabel); | |
1909 | grandmomPDG = grandmomP->GetPdgCode(); | |
1910 | ||
1911 | if(grandmomPDG==pdg) break; | |
1912 | ||
1913 | grandmomLabel = grandmomP->GetMother(); | |
1914 | ||
1915 | } | |
1916 | ||
1917 | if(grandmomPDG==pdg && grandmomP->GetNDaughters()==2) | |
1918 | { | |
1919 | AliAODMCParticle * d1 = (AliAODMCParticle *) mcparticles->At(grandmomP->GetDaughter(0)); | |
1920 | AliAODMCParticle * d2 = (AliAODMCParticle *) mcparticles->At(grandmomP->GetDaughter(1)); | |
1921 | if(d1->GetPdgCode() == 22 && d1->GetPdgCode() == 22) | |
1922 | { | |
1923 | asym = (d1->E()-d2->E())/grandmomP->E(); | |
0d360687 | 1924 | fDaughMom .SetPxPyPzE(d1->Px(),d1->Py(),d1->Pz(),d1->E()); |
1925 | fDaughMom2.SetPxPyPzE(d2->Px(),d2->Py(),d2->Pz(),d2->E()); | |
1926 | angle = fDaughMom.Angle(fDaughMom2.Vect()); | |
8e81c2cf | 1927 | } |
1928 | } | |
1929 | else | |
1930 | { | |
1931 | ok=kFALSE; | |
64be16b4 | 1932 | printf("AliMCAnalysisUtils::GetMCDecayAsymmetryForPDG(AOD) - mother with PDG %d, not found! \n",pdg); |
8e81c2cf | 1933 | } |
1934 | ||
1935 | } // good label | |
1936 | } | |
1937 | ||
1938 | ok = kTRUE; | |
1939 | ||
8e81c2cf | 1940 | } |
1941 | ||
8a2dbbff | 1942 | //_________________________________________________________________________________________________ |
1943 | Int_t AliMCAnalysisUtils::GetNDaughters(Int_t label, const AliCaloTrackReader* reader, Bool_t & ok) | |
64be16b4 | 1944 | { |
0d360687 | 1945 | // Return the the number of daughters of a given MC particle. |
64be16b4 | 1946 | |
1947 | if(reader->ReadStack()) | |
1948 | { | |
1949 | if(!reader->GetStack()) | |
1950 | { | |
1951 | if (fDebug >=0) | |
1952 | printf("AliMCAnalysisUtils::GetNDaughters() - Stack is not available, check analysis settings in configuration file, STOP!!\n"); | |
1953 | ||
1954 | ok=kFALSE; | |
1955 | return -1; | |
1956 | } | |
1957 | if(label >= 0 && label < reader->GetStack()->GetNtrack()) | |
1958 | { | |
1959 | TParticle * momP = reader->GetStack()->Particle(label); | |
1960 | ok=kTRUE; | |
1961 | return momP->GetNDaughters(); | |
1962 | } | |
1963 | else | |
1964 | { | |
1965 | ok = kFALSE; | |
1966 | return -1; | |
1967 | } | |
1968 | } | |
1969 | else if(reader->ReadAODMCParticles()) | |
1970 | { | |
1971 | TClonesArray* mcparticles = reader->GetAODMCParticles(); | |
1972 | if(!mcparticles) | |
1973 | { | |
1974 | if(fDebug >= 0) | |
1975 | printf("AliMCAnalysisUtils::GetNDaughters() - AODMCParticles is not available, check analysis settings in configuration file!!\n"); | |
1976 | ||
1977 | ok=kFALSE; | |
1978 | return -1; | |
1979 | } | |
1980 | ||
1981 | Int_t nprimaries = mcparticles->GetEntriesFast(); | |
1982 | if(label >= 0 && label < nprimaries) | |
1983 | { | |
1984 | AliAODMCParticle * momP = (AliAODMCParticle *) mcparticles->At(label); | |
1985 | ok = kTRUE; | |
1986 | return momP->GetNDaughters(); | |
1987 | } | |
1988 | else | |
1989 | { | |
1990 | ok = kFALSE; | |
1991 | return -1; | |
1992 | } | |
1993 | } | |
1994 | ||
1995 | ok = kFALSE; | |
1996 | ||
1997 | return -1; | |
1998 | } | |
1999 | ||
8a2dbbff | 2000 | //_________________________________________________________________________________ |
2001 | Int_t AliMCAnalysisUtils::GetNOverlaps(const Int_t * label, UInt_t nlabels, | |
2002 | Int_t mctag, Int_t mesonLabel, | |
4914e781 | 2003 | AliCaloTrackReader * reader, Int_t *overpdg) |
2004 | { | |
2005 | // Compare the primary depositing more energy with the rest, | |
0d360687 | 2006 | // if no photon/electron (conversion) or neutral meson as comon ancestor, consider it as other particle contributing. |
2007 | // Give as input the meson label in case it was a pi0 or eta merged cluster. | |
2008 | // Init overpdg with nlabels. | |
4914e781 | 2009 | |
2010 | Int_t ancPDG = 0, ancStatus = -1; | |
0d360687 | 2011 | TVector3 prodVertex; |
4914e781 | 2012 | Int_t ancLabel = 0; |
2013 | Int_t noverlaps = 0; | |
2014 | Bool_t ok = kFALSE; | |
2015 | ||
2016 | for (UInt_t ilab = 1; ilab < nlabels; ilab++ ) | |
2017 | { | |
0d360687 | 2018 | ancLabel = CheckCommonAncestor(label[0],label[ilab],reader,ancPDG,ancStatus,fMotherMom,prodVertex); |
4914e781 | 2019 | |
2020 | //printf("Overlaps, i %d: Main Label %d, second label %d, ancestor: Label %d, pdg %d - tag %d \n", | |
2021 | // ilab,label[0],label[ilab],ancLabel,ancPDG, mctag); | |
2022 | ||
2023 | Bool_t overlap = kFALSE; | |
2024 | ||
2025 | if ( ancLabel < 0 ) | |
2026 | { | |
2027 | overlap = kTRUE; | |
2028 | //printf("\t \t \t No Label = %d\n",ancLabel); | |
2029 | } | |
2030 | else if( ( ancPDG==111 || ancPDG==221 ) && ( CheckTagBit(mctag,kMCPi0) || CheckTagBit(mctag,kMCEta)) && mesonLabel != ancLabel) | |
2031 | { | |
2032 | //printf("\t \t meson Label %d, ancestor Label %d\n",mesonLabel,ancLabel); | |
2033 | overlap = kTRUE; | |
2034 | } | |
2035 | else if( ancPDG!=22 && TMath::Abs(ancPDG)!=11 && ancPDG != 111 && ancPDG != 221 ) | |
2036 | { | |
2037 | //printf("\t \t \t Non EM PDG = %d\n",ancPDG); | |
2038 | overlap = kTRUE ; | |
2039 | } | |
2040 | ||
2041 | if( !overlap ) continue ; | |
2042 | ||
2043 | // We have at least one overlap | |
2044 | ||
2045 | //printf("Overlap!!!!!!!!!!!!!!\n"); | |
2046 | ||
2047 | noverlaps++; | |
2048 | ||
2049 | // What is the origin of the overlap? | |
2050 | Bool_t mOK = 0, gOK = 0; | |
2051 | Int_t mpdg = -999999, gpdg = -1; | |
2052 | Int_t mstatus = -1, gstatus = -1; | |
2053 | Int_t gLabel = -1, ggLabel = -1; | |
0d360687 | 2054 | |
2055 | GetMother (label[ilab],reader,mpdg,mstatus,mOK); | |
2056 | fGMotherMom = | |
2057 | GetGrandMother(label[ilab],reader,gpdg,gstatus,gOK, gLabel,ggLabel); | |
4914e781 | 2058 | |
2059 | //printf("\t Overlap!, mother pdg %d; grand mother pdg %d",mpdg,gpdg); | |
2060 | ||
2061 | if( ( mpdg == 22 || TMath::Abs(mpdg==11) ) && | |
2062 | ( gpdg == 22 || TMath::Abs(gpdg==11) ) && | |
2063 | gLabel >=0 ) | |
2064 | { | |
2065 | Int_t labeltmp = gLabel; | |
2066 | while( ( gpdg == 22 || TMath::Abs(gpdg==11) ) && gLabel >=0 ) | |
2067 | { | |
2068 | mpdg=gpdg; | |
0d360687 | 2069 | fGMotherMom = GetGrandMother(labeltmp,reader,gpdg,gstatus,ok, gLabel,ggLabel); |
4914e781 | 2070 | labeltmp=gLabel; |
2071 | } | |
2072 | } | |
2073 | overpdg[noverlaps-1] = mpdg; | |
2074 | } | |
2075 | ||
2076 | return noverlaps ; | |
2077 | ||
2078 | } | |
2079 | ||
c5693f62 | 2080 | //________________________________________________________ |
7cd4e982 | 2081 | void AliMCAnalysisUtils::Print(const Option_t * opt) const |
2082 | { | |
0d360687 | 2083 | // Print some relevant parameters set for the analysis. |
f8006433 | 2084 | |
2085 | if(! opt) | |
2086 | return; | |
2087 | ||
2088 | printf("***** Print: %s %s ******\n", GetName(), GetTitle() ) ; | |
2089 | ||
2090 | printf("Debug level = %d\n",fDebug); | |
0d360687 | 2091 | printf("MC Generator = %s\n",fMCGeneratorString.Data()); |
f8006433 | 2092 | printf(" \n"); |
2093 | ||
7cd4e982 | 2094 | } |
2095 | ||
8a2dbbff | 2096 | //__________________________________________________ |
2097 | void AliMCAnalysisUtils::PrintMCTag(Int_t tag) const | |
bbb79e66 | 2098 | { |
0d360687 | 2099 | // Print the assigned origins to this particle. |
bbb79e66 | 2100 | |
4250012f | 2101 | printf("AliMCAnalysisUtils::PrintMCTag() - tag %d \n photon %d, conv %d, prompt %d, frag %d, isr %d, \n pi0 decay %d, eta decay %d, other decay %d pi0 %d, eta %d \n electron %d, muon %d,pion %d, proton %d, neutron %d, \n kaon %d, a-proton %d, a-neutron %d, unk %d, bad %d\n", |
bbb79e66 | 2102 | tag, |
2103 | CheckTagBit(tag,kMCPhoton), | |
bfc290d0 | 2104 | CheckTagBit(tag,kMCConversion), |
bbb79e66 | 2105 | CheckTagBit(tag,kMCPrompt), |
2106 | CheckTagBit(tag,kMCFragmentation), | |
2107 | CheckTagBit(tag,kMCISR), | |
2108 | CheckTagBit(tag,kMCPi0Decay), | |
2109 | CheckTagBit(tag,kMCEtaDecay), | |
2110 | CheckTagBit(tag,kMCOtherDecay), | |
bbb79e66 | 2111 | CheckTagBit(tag,kMCPi0), |
bfc290d0 | 2112 | CheckTagBit(tag,kMCEta), |
bbb79e66 | 2113 | CheckTagBit(tag,kMCElectron), |
bbb79e66 | 2114 | CheckTagBit(tag,kMCMuon), |
2115 | CheckTagBit(tag,kMCPion), | |
2116 | CheckTagBit(tag,kMCProton), | |
2117 | CheckTagBit(tag,kMCAntiNeutron), | |
2118 | CheckTagBit(tag,kMCKaon), | |
2119 | CheckTagBit(tag,kMCAntiProton), | |
2120 | CheckTagBit(tag,kMCAntiNeutron), | |
bfc290d0 | 2121 | CheckTagBit(tag,kMCUnknown), |
bbb79e66 | 2122 | CheckTagBit(tag,kMCBadLabel) |
2123 | ); | |
bbb79e66 | 2124 | } |
2125 | ||
0d360687 | 2126 | //__________________________________________________ |
2127 | void AliMCAnalysisUtils::SetMCGenerator(Int_t mcgen) | |
2128 | { | |
2129 | // Set the generator type. | |
2130 | ||
2131 | fMCGenerator = mcgen ; | |
2132 | if (mcgen == kPythia) fMCGeneratorString = "PYTHIA"; | |
2133 | else if(mcgen == kHerwig) fMCGeneratorString = "HERWIG"; | |
2134 | else if(mcgen == kHijing) fMCGeneratorString = "HIJING"; | |
2135 | else | |
2136 | { | |
2137 | fMCGeneratorString = ""; | |
2138 | fMCGenerator = kBoxLike ; | |
2139 | } | |
2140 | ||
2141 | } | |
2142 | ||
2143 | //__________________________________________________ | |
2144 | void AliMCAnalysisUtils::SetMCGenerator(TString mcgen) | |
2145 | { | |
2146 | // Set the generator type. | |
2147 | ||
2148 | fMCGeneratorString = mcgen ; | |
2149 | ||
2150 | if (mcgen == "PYTHIA") fMCGenerator = kPythia; | |
2151 | else if(mcgen == "HERWIG") fMCGenerator = kHerwig; | |
2152 | else if(mcgen == "HIJING") fMCGenerator = kHijing; | |
2153 | else | |
2154 | { | |
2155 | fMCGenerator = kBoxLike; | |
2156 | fMCGeneratorString = "" ; | |
2157 | } | |
2158 | } | |
2159 | ||
2160 | ||
7cd4e982 | 2161 |