]>
Commit | Line | Data |
---|---|---|
563113d0 | 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 | **************************************************************************/ | |
15 | ||
16 | //////////////////////////////////////////////////////////////////////////// | |
17 | // ---- CORRECTION FRAMEWORK ---- | |
18 | // class AliCFParticleGenCuts implementation | |
19 | // Using this class a user may define selections relative to | |
20 | // MC particle (AliMCParticle) using generation-level information. | |
21 | //////////////////////////////////////////////////////////////////////////// | |
22 | // author : R. Vernet (renaud.vernet@cern.ch) | |
23 | //////////////////////////////////////////////////////////////////////////// | |
24 | ||
25 | #include "AliLog.h" | |
26 | #include "AliCFParticleGenCuts.h" | |
27 | #include "TParticle.h" | |
28 | #include "TParticlePDG.h" | |
563113d0 | 29 | #include "AliMCEvent.h" |
30 | #include "TObject.h" | |
31 | #include "AliStack.h" | |
107a3100 | 32 | #include "TH1F.h" |
33 | #include "TH2F.h" | |
34 | #include "TBits.h" | |
35 | #include "TList.h" | |
c7803356 | 36 | #include "TArrayF.h" |
8fe1a43d | 37 | #include "TDecayChannel.h" |
b95f6a36 | 38 | #include "AliAODMCParticle.h" |
39 | #include "AliAODEvent.h" | |
563113d0 | 40 | |
41 | ClassImp(AliCFParticleGenCuts) | |
42 | ||
43 | //______________________________ | |
44 | AliCFParticleGenCuts::AliCFParticleGenCuts() : | |
45 | AliCFCutBase(), | |
b95f6a36 | 46 | fIsAODMC(0), |
563113d0 | 47 | fMCInfo(0x0), |
48 | fRequireIsCharged(0), | |
107a3100 | 49 | fRequireIsNeutral(0), |
563113d0 | 50 | fRequireIsPrimary(0), |
51 | fRequireIsSecondary(0), | |
52 | fRequirePdgCode(0), | |
3ff4a092 | 53 | fRequireAbsolutePdg(0), |
563113d0 | 54 | fPdgCode(0), |
55 | fProdVtxXMin (-1.e+09), | |
56 | fProdVtxYMin (-1.e+09), | |
57 | fProdVtxZMin (-1.e+09), | |
58 | fProdVtxXMax ( 1.e+09), | |
59 | fProdVtxYMax ( 1.e+09), | |
60 | fProdVtxZMax ( 1.e+09), | |
61 | fDecayVtxXMin(-1.e+09), | |
62 | fDecayVtxYMin(-1.e+09), | |
63 | fDecayVtxZMin(-1.e+09), | |
64 | fDecayVtxXMax( 1.e+09), | |
65 | fDecayVtxYMax( 1.e+09), | |
66 | fDecayVtxZMax( 1.e+09), | |
107a3100 | 67 | fDecayLengthMin(-1.), |
563113d0 | 68 | fDecayLengthMax(1.e+09), |
107a3100 | 69 | fDecayRxyMin(-1), |
70 | fDecayRxyMax(1.e+09), | |
8fe1a43d | 71 | fDecayChannel(0x0), |
107a3100 | 72 | fhCutStatistics(0x0), |
73 | fhCutCorrelation(0x0), | |
c7803356 | 74 | fCutValues(new TArrayF(kNCuts)), |
107a3100 | 75 | fBitmap(new TBits(0)) |
563113d0 | 76 | { |
77 | // | |
78 | //ctor | |
79 | // | |
107a3100 | 80 | for (int i=0; i<kNCuts; i++) |
81 | for (int j=0; j<kNStepQA; j++) | |
82 | fhQA[i][j]=0x0; | |
563113d0 | 83 | } |
84 | ||
85 | //______________________________ | |
86 | AliCFParticleGenCuts::AliCFParticleGenCuts(const Char_t* name, const Char_t* title) : | |
87 | AliCFCutBase(name,title), | |
b95f6a36 | 88 | fIsAODMC(0), |
563113d0 | 89 | fMCInfo(0x0), |
90 | fRequireIsCharged(0), | |
107a3100 | 91 | fRequireIsNeutral(0), |
563113d0 | 92 | fRequireIsPrimary(0), |
93 | fRequireIsSecondary(0), | |
94 | fRequirePdgCode(0), | |
3ff4a092 | 95 | fRequireAbsolutePdg(0), |
563113d0 | 96 | fPdgCode(0), |
97 | fProdVtxXMin (-1.e+09), | |
98 | fProdVtxYMin (-1.e+09), | |
99 | fProdVtxZMin (-1.e+09), | |
100 | fProdVtxXMax ( 1.e+09), | |
101 | fProdVtxYMax ( 1.e+09), | |
102 | fProdVtxZMax ( 1.e+09), | |
103 | fDecayVtxXMin(-1.e+09), | |
104 | fDecayVtxYMin(-1.e+09), | |
105 | fDecayVtxZMin(-1.e+09), | |
106 | fDecayVtxXMax( 1.e+09), | |
107 | fDecayVtxYMax( 1.e+09), | |
108 | fDecayVtxZMax( 1.e+09), | |
107a3100 | 109 | fDecayLengthMin(-1.), |
563113d0 | 110 | fDecayLengthMax(1.e+09), |
107a3100 | 111 | fDecayRxyMin(-1.), |
112 | fDecayRxyMax(1.e+09), | |
8fe1a43d | 113 | fDecayChannel(0x0), |
107a3100 | 114 | fhCutStatistics(0x0), |
115 | fhCutCorrelation(0x0), | |
c7803356 | 116 | fCutValues(new TArrayF(kNCuts)), |
107a3100 | 117 | fBitmap(new TBits(0)) |
563113d0 | 118 | { |
119 | // | |
120 | //ctor | |
121 | // | |
107a3100 | 122 | for (int i=0; i<kNCuts; i++) |
123 | for (int j=0; j<kNStepQA; j++) | |
124 | fhQA[i][j]=0x0; | |
563113d0 | 125 | } |
126 | ||
127 | //______________________________ | |
128 | AliCFParticleGenCuts::AliCFParticleGenCuts(const AliCFParticleGenCuts& c) : | |
129 | AliCFCutBase(c), | |
b95f6a36 | 130 | fIsAODMC(c.fIsAODMC), |
563113d0 | 131 | fMCInfo(c.fMCInfo), |
132 | fRequireIsCharged(c.fRequireIsCharged), | |
107a3100 | 133 | fRequireIsNeutral(c.fRequireIsNeutral), |
563113d0 | 134 | fRequireIsPrimary(c.fRequireIsPrimary), |
135 | fRequireIsSecondary(c.fRequireIsSecondary), | |
136 | fRequirePdgCode(c.fRequirePdgCode), | |
3ff4a092 | 137 | fRequireAbsolutePdg(c.fRequireAbsolutePdg), |
563113d0 | 138 | fPdgCode(c.fPdgCode), |
139 | fProdVtxXMin (c.fProdVtxXMin), | |
140 | fProdVtxYMin (c.fProdVtxYMin), | |
141 | fProdVtxZMin (c.fProdVtxZMin), | |
142 | fProdVtxXMax (c.fProdVtxXMax), | |
143 | fProdVtxYMax (c.fProdVtxYMax), | |
144 | fProdVtxZMax (c.fProdVtxZMax), | |
145 | fDecayVtxXMin(c.fDecayVtxXMin), | |
146 | fDecayVtxYMin(c.fDecayVtxYMin), | |
147 | fDecayVtxZMin(c.fDecayVtxZMin), | |
148 | fDecayVtxXMax(c.fDecayVtxXMax), | |
149 | fDecayVtxYMax(c.fDecayVtxYMax), | |
150 | fDecayVtxZMax(c.fDecayVtxZMax), | |
151 | fDecayLengthMin(c.fDecayLengthMin), | |
152 | fDecayLengthMax(c.fDecayLengthMin), | |
153 | fDecayRxyMin(c.fDecayLengthMin), | |
107a3100 | 154 | fDecayRxyMax(c.fDecayLengthMin), |
8fe1a43d | 155 | fDecayChannel(c.fDecayChannel), |
107a3100 | 156 | fhCutStatistics(new TH1F(*c.fhCutStatistics)), |
157 | fhCutCorrelation(new TH2F(*c.fhCutCorrelation)), | |
c7803356 | 158 | fCutValues(new TArrayF(*c.fCutValues)), |
107a3100 | 159 | fBitmap(new TBits(*c.fBitmap)) |
563113d0 | 160 | { |
161 | // | |
162 | //copy ctor | |
163 | // | |
107a3100 | 164 | for (int i=0; i<kNCuts; i++) |
165 | for (int j=0; j<kNStepQA; j++) | |
166 | fhQA[i][j]=(TH1F*)c.fhQA[i][j]->Clone(); | |
563113d0 | 167 | } |
168 | ||
169 | //______________________________ | |
170 | AliCFParticleGenCuts& AliCFParticleGenCuts::operator=(const AliCFParticleGenCuts& c) | |
171 | { | |
172 | // | |
173 | // Assignment operator | |
174 | // | |
175 | if (this != &c) { | |
176 | AliCFCutBase::operator=(c) ; | |
b95f6a36 | 177 | fIsAODMC=c.fIsAODMC; |
563113d0 | 178 | fMCInfo=c.fMCInfo; |
179 | fRequireIsCharged=c.fRequireIsCharged; | |
107a3100 | 180 | fRequireIsNeutral=c.fRequireIsNeutral; |
563113d0 | 181 | fRequireIsPrimary=c.fRequireIsPrimary; |
182 | fRequireIsSecondary=c.fRequireIsSecondary; | |
183 | fRequirePdgCode=c.fRequirePdgCode; | |
3ff4a092 | 184 | fRequireAbsolutePdg=c.fRequireAbsolutePdg; |
563113d0 | 185 | fPdgCode=c.fPdgCode; |
186 | fProdVtxXMin=c.fProdVtxXMin; | |
187 | fProdVtxYMin=c.fProdVtxYMin; | |
188 | fProdVtxZMin=c.fProdVtxZMin; | |
189 | fProdVtxXMax=c.fProdVtxXMax; | |
190 | fProdVtxYMax=c.fProdVtxYMax; | |
191 | fProdVtxZMax=c.fProdVtxZMax; | |
192 | fDecayVtxXMin=c.fDecayVtxXMin; | |
193 | fDecayVtxYMin=c.fDecayVtxYMin; | |
194 | fDecayVtxZMin=c.fDecayVtxZMin; | |
195 | fDecayVtxXMax=c.fDecayVtxXMax; | |
196 | fDecayVtxYMax=c.fDecayVtxYMax; | |
197 | fDecayVtxZMax=c.fDecayVtxZMax; | |
198 | fDecayLengthMin=c.fDecayVtxZMax; | |
199 | fDecayLengthMax=c.fDecayLengthMax; | |
200 | fDecayRxyMin=c.fDecayRxyMin; | |
201 | fDecayRxyMax=c.fDecayRxyMax; | |
8fe1a43d | 202 | fDecayChannel=c.fDecayChannel; |
c7803356 | 203 | fCutValues=new TArrayF(*c.fCutValues); |
107a3100 | 204 | fBitmap=new TBits(*c.fBitmap); |
205 | ||
206 | if (fhCutStatistics) fhCutStatistics =new TH1F(*c.fhCutStatistics) ; | |
207 | if (fhCutCorrelation) fhCutCorrelation=new TH2F(*c.fhCutCorrelation); | |
208 | ||
209 | for (int i=0; i<kNCuts; i++) | |
210 | for (int j=0; j<kNStepQA; j++) | |
211 | fhQA[i][j]=(TH1F*)c.fhQA[i][j]->Clone(); | |
563113d0 | 212 | } |
213 | return *this ; | |
214 | } | |
215 | ||
216 | //______________________________ | |
217 | Bool_t AliCFParticleGenCuts::IsSelected(TObject* obj) { | |
218 | // | |
219 | // check if selections on 'obj' are passed | |
220 | // 'obj' must be an AliMCParticle | |
221 | // | |
222 | ||
b95f6a36 | 223 | if (!obj) return kFALSE ; |
224 | ||
225 | if (!fIsAODMC) SelectionBitMap((AliMCParticle*) obj); | |
226 | else SelectionBitMap((AliAODMCParticle*)obj); | |
107a3100 | 227 | |
228 | if (fIsQAOn) FillHistograms(obj,0); | |
229 | ||
230 | for (UInt_t icut=0; icut<fBitmap->GetNbits();icut++) | |
231 | if (!fBitmap->TestBitNumber(icut)) return kFALSE ; | |
232 | ||
233 | if (fIsQAOn) FillHistograms(obj,1); | |
234 | return kTRUE; | |
235 | } | |
c7803356 | 236 | |
107a3100 | 237 | //__________________________________________________________________________________ |
b95f6a36 | 238 | void AliCFParticleGenCuts::SelectionBitMap(AliMCParticle* mcPart) |
107a3100 | 239 | { |
240 | // | |
241 | // test if the track passes the single cuts | |
242 | // and store the information in a bitmap | |
243 | // | |
b95f6a36 | 244 | |
c7803356 | 245 | for (UInt_t i=0; i<kNCuts; i++) { |
246 | fBitmap->SetBitNumber(i,kFALSE); | |
247 | fCutValues->SetAt((Double32_t)0,i) ; | |
248 | } | |
107a3100 | 249 | |
c7803356 | 250 | // fill the cut array |
b95f6a36 | 251 | Double32_t partVx=(Double32_t)mcPart->Xv(); |
252 | Double32_t partVy=(Double32_t)mcPart->Yv(); | |
253 | Double32_t partVz=(Double32_t)mcPart->Zv(); | |
107a3100 | 254 | |
107a3100 | 255 | Double32_t decayVx=0.; |
256 | Double32_t decayVy=0.; | |
257 | Double32_t decayVz=0.; | |
258 | Double32_t decayL=0.; | |
259 | Double32_t decayRxy=0.; | |
260 | ||
b95f6a36 | 261 | TParticle* part = mcPart->Particle(); |
262 | AliStack* stack = ((AliMCEvent*)fMCInfo)->Stack(); | |
263 | TParticle* daughter=0x0; | |
107a3100 | 264 | if ( part->GetNDaughters() > 0 ) { |
265 | daughter = stack->Particle(part->GetFirstDaughter()) ; | |
266 | decayVx=(Double32_t)daughter->Vx(); | |
267 | decayVy=(Double32_t)daughter->Vy(); | |
268 | decayVz=(Double32_t)daughter->Vz(); | |
269 | decayL = TMath::Sqrt(TMath::Power(partVx-decayVx,2) + | |
270 | TMath::Power(partVy-decayVy,2) + | |
271 | TMath::Power(partVz-decayVz,2) ) ; | |
272 | decayRxy = TMath::Sqrt(TMath::Power(decayVx,2) + TMath::Power(decayVy,2) ) ; | |
273 | } | |
274 | ||
c7803356 | 275 | fCutValues->SetAt(partVx ,kCutProdVtxXMin) ; |
276 | fCutValues->SetAt(partVy ,kCutProdVtxYMin) ; | |
277 | fCutValues->SetAt(partVz ,kCutProdVtxZMin) ; | |
278 | fCutValues->SetAt(partVx ,kCutProdVtxXMax) ; | |
279 | fCutValues->SetAt(partVy ,kCutProdVtxYMax) ; | |
280 | fCutValues->SetAt(partVz ,kCutProdVtxZMax) ; | |
281 | fCutValues->SetAt(decayVx ,kCutDecVtxXMin) ; | |
282 | fCutValues->SetAt(decayVy ,kCutDecVtxYMin) ; | |
283 | fCutValues->SetAt(decayVz ,kCutDecVtxZMin) ; | |
284 | fCutValues->SetAt(decayVx ,kCutDecVtxXMax) ; | |
285 | fCutValues->SetAt(decayVy ,kCutDecVtxYMax) ; | |
286 | fCutValues->SetAt(decayVz ,kCutDecVtxZMax) ; | |
287 | fCutValues->SetAt(decayL ,kCutDecLgthMin) ; | |
288 | fCutValues->SetAt(decayL ,kCutDecLgthMax) ; | |
289 | fCutValues->SetAt(decayRxy,kCutDecRxyMin) ; | |
290 | fCutValues->SetAt(decayRxy,kCutDecRxyMax) ; | |
291 | ||
107a3100 | 292 | // cut on charge |
293 | if ( fRequireIsCharged || fRequireIsNeutral ) { | |
c7803356 | 294 | if (fRequireIsCharged && IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ; |
295 | if (fRequireIsNeutral && !IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ; | |
563113d0 | 296 | } |
c7803356 | 297 | else fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ; |
563113d0 | 298 | |
107a3100 | 299 | // cut on primary/secondary |
300 | if ( fRequireIsPrimary || fRequireIsSecondary) { | |
b95f6a36 | 301 | if (fRequireIsPrimary && IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec); |
302 | if (fRequireIsSecondary && !IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec); | |
563113d0 | 303 | } |
c7803356 | 304 | else fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec); |
305 | ||
107a3100 | 306 | // cut on PDG code |
c7803356 | 307 | if ( fRequirePdgCode ) { |
3ff4a092 | 308 | if (IsA(mcPart,fPdgCode,fRequireAbsolutePdg)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode); |
563113d0 | 309 | } |
c7803356 | 310 | else fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode); |
107a3100 | 311 | |
8fe1a43d | 312 | // cut on decay channel |
313 | if ( fDecayChannel ) { | |
314 | Bool_t goodDecay = kTRUE ; | |
315 | Short_t nDaughters = mcPart->Particle()->GetNDaughters() ; | |
316 | if (nDaughters != fDecayChannel->NDaughters()) goodDecay = kFALSE ; | |
2f8e8cad | 317 | //now number of daughters is OK |
8fe1a43d | 318 | if (goodDecay) { |
319 | // now check if decay channel is respected | |
320 | // first try | |
321 | for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) { | |
322 | TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(iDaughter)) ; | |
323 | if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;} | |
324 | } | |
325 | if (!goodDecay) { | |
326 | //second try inverting the order of the daughters | |
327 | goodDecay = kTRUE ; | |
328 | for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) { | |
329 | TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(nDaughters-(iDaughter+1))) ; | |
330 | if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;} | |
331 | } | |
332 | } | |
2f8e8cad | 333 | if (!goodDecay && fRequireAbsolutePdg) { |
334 | //now tries inverting the sign of the daughters in case the anti-particle is also looked at | |
335 | // third try | |
336 | goodDecay = kTRUE ; | |
337 | for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) { | |
338 | TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(iDaughter)) ; | |
339 | if (daug->GetPdgCode() != -fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;} | |
340 | } | |
341 | if (!goodDecay) { | |
342 | //fourth try inverting the order of the daughters | |
343 | goodDecay = kTRUE ; | |
344 | for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) { | |
345 | TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(nDaughters-(iDaughter+1))) ; | |
346 | if (daug->GetPdgCode() != -fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;} | |
347 | } | |
348 | } | |
349 | } //end check anti-particle | |
350 | } //end # daughters OK | |
8fe1a43d | 351 | fCutValues->SetAt((Double32_t)goodDecay,kCutDecayChannel) ; |
2f8e8cad | 352 | } //end require decay channel |
8fe1a43d | 353 | else fCutValues->SetAt((Double32_t)kTRUE,kCutDecayChannel); |
354 | ||
107a3100 | 355 | |
c7803356 | 356 | // now array of cut is build, fill the bitmap consequently |
357 | Int_t iCutBit = -1; | |
358 | if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
359 | if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
360 | if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
361 | if ( fCutValues->At(++iCutBit) > fProdVtxXMin) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
362 | if ( fCutValues->At(++iCutBit) < fProdVtxXMax) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
363 | if ( fCutValues->At(++iCutBit) > fProdVtxYMin) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
364 | if ( fCutValues->At(++iCutBit) < fProdVtxYMax) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
365 | if ( fCutValues->At(++iCutBit) > fProdVtxZMin) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
366 | if ( fCutValues->At(++iCutBit) < fProdVtxZMax) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
367 | if ( fCutValues->At(++iCutBit) > fDecayVtxXMin) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
368 | if ( fCutValues->At(++iCutBit) < fDecayVtxXMax) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
369 | if ( fCutValues->At(++iCutBit) > fDecayVtxYMin) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
370 | if ( fCutValues->At(++iCutBit) < fDecayVtxYMax) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
371 | if ( fCutValues->At(++iCutBit) > fDecayVtxZMin) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
372 | if ( fCutValues->At(++iCutBit) < fDecayVtxZMax) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
373 | if ( fCutValues->At(++iCutBit) > fDecayLengthMin) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
374 | if ( fCutValues->At(++iCutBit) < fDecayLengthMax) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
375 | if ( fCutValues->At(++iCutBit) > fDecayRxyMin) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
376 | if ( fCutValues->At(++iCutBit) < fDecayRxyMax) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
8fe1a43d | 377 | if ( fCutValues->At(++iCutBit) != 0 ) fBitmap->SetBitNumber(iCutBit,kTRUE); |
107a3100 | 378 | } |
379 | ||
b95f6a36 | 380 | //__________________________________________________________________________________ |
381 | void AliCFParticleGenCuts::SelectionBitMap(AliAODMCParticle* mcPart) | |
382 | { | |
383 | // | |
384 | // test if the track passes the single cuts | |
385 | // and store the information in a bitmap | |
386 | // | |
387 | ||
388 | for (UInt_t i=0; i<kNCuts; i++) { | |
389 | fBitmap->SetBitNumber(i,kFALSE); | |
390 | fCutValues->SetAt((Double32_t)0,i) ; | |
391 | } | |
392 | ||
393 | // fill the cut array | |
394 | Double32_t partVx=(Double32_t)mcPart->Xv(); | |
395 | Double32_t partVy=(Double32_t)mcPart->Yv(); | |
396 | Double32_t partVz=(Double32_t)mcPart->Zv(); | |
397 | ||
398 | Double32_t decayVx=0.; | |
399 | Double32_t decayVy=0.; | |
400 | Double32_t decayVz=0.; | |
401 | Double32_t decayL=0.; | |
402 | Double32_t decayRxy=0.; | |
403 | ||
404 | AliAODEvent* aod = dynamic_cast<AliAODEvent*>(fMCInfo); | |
405 | TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName())); | |
406 | AliAODMCParticle* daughter = 0x0 ; | |
407 | ||
408 | if (mcPart->GetDaughter(0)>0) { | |
409 | daughter = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(0))); | |
410 | ||
411 | decayVx=(Double32_t)daughter->Xv(); | |
412 | decayVy=(Double32_t)daughter->Yv(); | |
413 | decayVz=(Double32_t)daughter->Zv(); | |
414 | decayL = TMath::Sqrt(TMath::Power(partVx-decayVx,2) + | |
415 | TMath::Power(partVy-decayVy,2) + | |
416 | TMath::Power(partVz-decayVz,2) ) ; | |
417 | decayRxy = TMath::Sqrt(TMath::Power(decayVx,2) + TMath::Power(decayVy,2) ) ; | |
418 | ||
419 | } | |
420 | ||
421 | fCutValues->SetAt(partVx ,kCutProdVtxXMin) ; | |
422 | fCutValues->SetAt(partVy ,kCutProdVtxYMin) ; | |
423 | fCutValues->SetAt(partVz ,kCutProdVtxZMin) ; | |
424 | fCutValues->SetAt(partVx ,kCutProdVtxXMax) ; | |
425 | fCutValues->SetAt(partVy ,kCutProdVtxYMax) ; | |
426 | fCutValues->SetAt(partVz ,kCutProdVtxZMax) ; | |
427 | fCutValues->SetAt(decayVx ,kCutDecVtxXMin) ; | |
428 | fCutValues->SetAt(decayVy ,kCutDecVtxYMin) ; | |
429 | fCutValues->SetAt(decayVz ,kCutDecVtxZMin) ; | |
430 | fCutValues->SetAt(decayVx ,kCutDecVtxXMax) ; | |
431 | fCutValues->SetAt(decayVy ,kCutDecVtxYMax) ; | |
432 | fCutValues->SetAt(decayVz ,kCutDecVtxZMax) ; | |
433 | fCutValues->SetAt(decayL ,kCutDecLgthMin) ; | |
434 | fCutValues->SetAt(decayL ,kCutDecLgthMax) ; | |
435 | fCutValues->SetAt(decayRxy,kCutDecRxyMin) ; | |
436 | fCutValues->SetAt(decayRxy,kCutDecRxyMax) ; | |
437 | ||
438 | // cut on charge | |
439 | if ( fRequireIsCharged || fRequireIsNeutral ) { | |
440 | if (fRequireIsCharged && IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ; | |
441 | if (fRequireIsNeutral && !IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ; | |
442 | } | |
443 | else fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ; | |
444 | ||
445 | // cut on primary/secondary | |
446 | if ( fRequireIsPrimary || fRequireIsSecondary) { | |
447 | if (fRequireIsPrimary && IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec); | |
448 | if (fRequireIsSecondary && !IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec); | |
449 | } | |
450 | else fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec); | |
451 | ||
452 | // cut on PDG code | |
453 | if ( fRequirePdgCode ) { | |
3ff4a092 | 454 | if (IsA(mcPart,fPdgCode,fRequireAbsolutePdg)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode); |
b95f6a36 | 455 | } |
456 | else fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode); | |
457 | ||
458 | // cut on decay channel | |
459 | if ( fDecayChannel ) { | |
460 | Bool_t goodDecay = kTRUE ; | |
461 | Short_t nDaughters = 0 ; | |
462 | if (mcPart->GetDaughter(0) >=0) nDaughters = 1 + mcPart->GetDaughter(1) - mcPart->GetDaughter(0); | |
463 | ||
464 | if (nDaughters != fDecayChannel->NDaughters()) goodDecay = kFALSE ; | |
465 | if (goodDecay) { | |
466 | // now check if decay channel is respected | |
467 | // first try | |
468 | for (Int_t iDaughter = 0 ; iDaughter<nDaughters; iDaughter++) { | |
469 | AliAODMCParticle* daug = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(0)+iDaughter)); | |
470 | if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;} | |
471 | } | |
472 | if (!goodDecay) { | |
473 | //second try inverting the order of the daughters | |
474 | goodDecay = kTRUE ; | |
475 | for (Int_t iDaughter = 0 ; iDaughter<nDaughters; iDaughter++) { | |
476 | AliAODMCParticle* daug = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(1)-iDaughter)); | |
477 | if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;} | |
478 | } | |
479 | } | |
480 | } | |
481 | fCutValues->SetAt((Double32_t)goodDecay,kCutDecayChannel) ; | |
482 | } | |
483 | else fCutValues->SetAt((Double32_t)kTRUE,kCutDecayChannel); | |
484 | ||
485 | ||
486 | // now array of cut is build, fill the bitmap consequently | |
487 | Int_t iCutBit = -1; | |
488 | if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
489 | if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
490 | if ( fCutValues->At(++iCutBit) !=0 ) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
491 | if ( fCutValues->At(++iCutBit) > fProdVtxXMin) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
492 | if ( fCutValues->At(++iCutBit) < fProdVtxXMax) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
493 | if ( fCutValues->At(++iCutBit) > fProdVtxYMin) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
494 | if ( fCutValues->At(++iCutBit) < fProdVtxYMax) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
495 | if ( fCutValues->At(++iCutBit) > fProdVtxZMin) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
496 | if ( fCutValues->At(++iCutBit) < fProdVtxZMax) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
497 | if ( fCutValues->At(++iCutBit) > fDecayVtxXMin) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
498 | if ( fCutValues->At(++iCutBit) < fDecayVtxXMax) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
499 | if ( fCutValues->At(++iCutBit) > fDecayVtxYMin) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
500 | if ( fCutValues->At(++iCutBit) < fDecayVtxYMax) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
501 | if ( fCutValues->At(++iCutBit) > fDecayVtxZMin) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
502 | if ( fCutValues->At(++iCutBit) < fDecayVtxZMax) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
503 | if ( fCutValues->At(++iCutBit) > fDecayLengthMin) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
504 | if ( fCutValues->At(++iCutBit) < fDecayLengthMax) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
505 | if ( fCutValues->At(++iCutBit) > fDecayRxyMin) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
506 | if ( fCutValues->At(++iCutBit) < fDecayRxyMax) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
507 | if ( fCutValues->At(++iCutBit) != 0 ) fBitmap->SetBitNumber(iCutBit,kTRUE); | |
508 | } | |
509 | ||
c7803356 | 510 | |
107a3100 | 511 | //__________________________________________________________________________________ |
033789c9 | 512 | void AliCFParticleGenCuts::FillHistograms(TObject* /*obj*/, Bool_t afterCuts) |
107a3100 | 513 | { |
514 | // | |
515 | // fill the QA histograms | |
516 | // | |
517 | ||
518 | for (int iCutNumber = 0; iCutNumber < kNCuts; iCutNumber++) | |
c7803356 | 519 | fhQA[iCutNumber][afterCuts]->Fill(fCutValues->At(iCutNumber)); |
107a3100 | 520 | |
521 | // fill cut statistics and cut correlation histograms with information from the bitmap | |
522 | if (afterCuts) return; | |
523 | ||
524 | // Number of single cuts in this class | |
525 | UInt_t ncuts = fBitmap->GetNbits(); | |
526 | for(UInt_t bit=0; bit<ncuts;bit++) { | |
c7803356 | 527 | if (!fBitmap->TestBitNumber(bit)) { |
107a3100 | 528 | fhCutStatistics->Fill(bit+1); |
529 | for (UInt_t bit2=bit; bit2<ncuts;bit2++) { | |
530 | if (!fBitmap->TestBitNumber(bit2)) | |
531 | fhCutCorrelation->Fill(bit+1,bit2+1); | |
532 | } | |
533 | } | |
563113d0 | 534 | } |
107a3100 | 535 | } |
536 | ||
537 | //__________________________________________________________________________________ | |
538 | void AliCFParticleGenCuts::AddQAHistograms(TList *qaList) { | |
539 | // | |
540 | // saves the histograms in a TList | |
541 | // | |
563113d0 | 542 | |
107a3100 | 543 | DefineHistograms(); |
563113d0 | 544 | |
107a3100 | 545 | qaList->Add(fhCutStatistics); |
546 | qaList->Add(fhCutCorrelation); | |
547 | ||
548 | for (Int_t j=0; j<kNStepQA; j++) { | |
549 | for(Int_t i=0; i<kNCuts; i++) | |
550 | qaList->Add(fhQA[i][j]); | |
551 | } | |
552 | } | |
553 | ||
554 | //__________________________________________________________________________________ | |
555 | void AliCFParticleGenCuts::DefineHistograms() { | |
556 | // | |
557 | // histograms for cut variables, cut statistics and cut correlations | |
558 | // | |
559 | Int_t color = 2; | |
560 | ||
561 | // book cut statistics and cut correlation histograms | |
562 | fhCutStatistics = new TH1F(Form("%s_cut_statistics",GetName()),"",kNCuts,0.5,kNCuts+0.5); | |
563 | fhCutStatistics->SetLineWidth(2); | |
564 | int k = 1; | |
565 | fhCutStatistics->GetXaxis()->SetBinLabel(k,"charge") ; k++; | |
566 | fhCutStatistics->GetXaxis()->SetBinLabel(k,"prim/sec") ; k++; | |
567 | fhCutStatistics->GetXaxis()->SetBinLabel(k,"PDG") ; k++; | |
568 | fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxXMin") ; k++; | |
569 | fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxXMax") ; k++; | |
570 | fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxYMin") ; k++; | |
571 | fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxYMax") ; k++; | |
572 | fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxZMin") ; k++; | |
573 | fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMax") ; k++; | |
574 | fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecXMin") ; k++; | |
575 | fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecXMax") ; k++; | |
576 | fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecYMin") ; k++; | |
577 | fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecYMax") ; k++; | |
578 | fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMin") ; k++; | |
579 | fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMax") ; k++; | |
580 | fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecLgthMin") ; k++; | |
581 | fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecLgthMax") ; k++; | |
582 | fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMin") ; k++; | |
583 | fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMax") ; k++; | |
8fe1a43d | 584 | fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecChannel") ; k++; |
107a3100 | 585 | |
586 | ||
587 | fhCutCorrelation = new TH2F(Form("%s_cut_correlation",GetName()),"",kNCuts,0.5,kNCuts+0.5,kNCuts,0.5,kNCuts+0.5); | |
588 | fhCutCorrelation->SetLineWidth(2); | |
589 | for (k=1; k<=kNCuts; k++) { | |
590 | fhCutCorrelation->GetXaxis()->SetBinLabel(k,fhCutStatistics->GetXaxis()->GetBinLabel(k)); | |
591 | fhCutCorrelation->GetYaxis()->SetBinLabel(k,fhCutStatistics->GetXaxis()->GetBinLabel(k)); | |
592 | } | |
593 | ||
594 | Char_t str[256]; | |
595 | for (int i=0; i<kNStepQA; i++) { | |
596 | if (i==0) sprintf(str," "); | |
597 | else sprintf(str,"_cut"); | |
8fe1a43d | 598 | fhQA[kCutCharge] [i] = new TH1F(Form("%s_charge%s" ,GetName(),str),"",2,0,2); |
599 | fhQA[kCutPrimSec] [i] = new TH1F(Form("%s_primSec%s" ,GetName(),str),"",2,0,2); | |
600 | fhQA[kCutPDGCode] [i] = new TH1F(Form("%s_pdgCode%s" ,GetName(),str),"",2,0,2); | |
601 | fhQA[kCutProdVtxXMin] [i] = new TH1F(Form("%s_prodVtxXMin%s" ,GetName(),str),"",100,0,10); | |
602 | fhQA[kCutProdVtxXMax] [i] = new TH1F(Form("%s_prodVtxXMax%s" ,GetName(),str),"",100,0,10); | |
603 | fhQA[kCutProdVtxYMin] [i] = new TH1F(Form("%s_prodVtxYMin%s" ,GetName(),str),"",100,0,10); | |
604 | fhQA[kCutProdVtxYMax] [i] = new TH1F(Form("%s_prodVtxYMax%s" ,GetName(),str),"",100,0,10); | |
605 | fhQA[kCutProdVtxZMin] [i] = new TH1F(Form("%s_prodVtxZMin%s" ,GetName(),str),"",100,0,10); | |
606 | fhQA[kCutProdVtxZMax] [i] = new TH1F(Form("%s_prodVtxZMax%s" ,GetName(),str),"",100,0,10); | |
607 | fhQA[kCutDecVtxXMin] [i] = new TH1F(Form("%s_decVtxXMin%s" ,GetName(),str),"",100,0,10); | |
608 | fhQA[kCutDecVtxXMax] [i] = new TH1F(Form("%s_decVtxXMax%s" ,GetName(),str),"",100,0,10); | |
609 | fhQA[kCutDecVtxYMin] [i] = new TH1F(Form("%s_decVtxYMin%s" ,GetName(),str),"",100,0,10); | |
610 | fhQA[kCutDecVtxYMax] [i] = new TH1F(Form("%s_decVtxYMax%s" ,GetName(),str),"",100,0,10); | |
611 | fhQA[kCutDecVtxZMin] [i] = new TH1F(Form("%s_decVtxZMin%s" ,GetName(),str),"",100,0,10); | |
612 | fhQA[kCutDecVtxZMax] [i] = new TH1F(Form("%s_decVtxZMax%s" ,GetName(),str),"",100,0,10); | |
613 | fhQA[kCutDecLgthMin] [i] = new TH1F(Form("%s_decLengthMin%s",GetName(),str),"",100,0,10); | |
614 | fhQA[kCutDecLgthMax] [i] = new TH1F(Form("%s_decLengthMax%s",GetName(),str),"",100,0,10); | |
615 | fhQA[kCutDecRxyMin] [i] = new TH1F(Form("%s_decRxyMin%s" ,GetName(),str),"",100,0,10); | |
616 | fhQA[kCutDecRxyMax] [i] = new TH1F(Form("%s_decRxyMax%s" ,GetName(),str),"",100,0,10); | |
617 | fhQA[kCutDecayChannel][i] = new TH1F(Form("%s_decayChannel%s",GetName(),str),"",2,0,2); | |
107a3100 | 618 | } |
619 | for(Int_t i=0; i<kNCuts; i++) fhQA[i][1]->SetLineColor(color); | |
563113d0 | 620 | } |
107a3100 | 621 | |
622 | ||
563113d0 | 623 | //______________________________ |
b95f6a36 | 624 | Bool_t AliCFParticleGenCuts::IsCharged(AliVParticle *mcPart) { |
563113d0 | 625 | // |
626 | //check if particle is charged. | |
627 | // | |
b95f6a36 | 628 | if (mcPart->Charge()==0) return kFALSE; |
563113d0 | 629 | return kTRUE; |
630 | } | |
631 | //______________________________ | |
b95f6a36 | 632 | Bool_t AliCFParticleGenCuts::IsPrimary(AliMCParticle *mcPart) { |
563113d0 | 633 | // |
634 | //check if particle is primary (standard definition) | |
635 | // | |
b95f6a36 | 636 | |
637 | AliStack* stack = ((AliMCEvent*)fMCInfo)->Stack(); | |
638 | ||
639 | if (!stack->IsPhysicalPrimary(mcPart->GetLabel())) return kFALSE; | |
563113d0 | 640 | return kTRUE; |
641 | } | |
642 | //______________________________ | |
b95f6a36 | 643 | Bool_t AliCFParticleGenCuts::IsPrimary(AliAODMCParticle *mcPart) { |
644 | // | |
645 | //check if particle is primary (standard definition) | |
646 | // | |
647 | ||
648 | if (!mcPart->IsPhysicalPrimary()) return kFALSE; | |
649 | return kTRUE; | |
650 | } | |
651 | //______________________________ | |
652 | Bool_t AliCFParticleGenCuts::IsPrimaryCharged(AliVParticle *mcPart) { | |
563113d0 | 653 | // |
654 | //check if a charged particle is primary (standard definition) | |
655 | // | |
b95f6a36 | 656 | |
657 | if (!fIsAODMC) { | |
658 | if (!IsPrimary((AliMCParticle*)mcPart) || !IsCharged(mcPart)) return kFALSE ; | |
659 | } | |
660 | else { | |
661 | if (!IsPrimary((AliAODMCParticle*)mcPart) || !IsCharged(mcPart)) return kFALSE ; | |
662 | } | |
563113d0 | 663 | return kTRUE; |
664 | } | |
665 | //______________________________ | |
666 | Bool_t AliCFParticleGenCuts::IsA(AliMCParticle *mcPart, Int_t pdg, Bool_t abs) { | |
667 | // | |
668 | //Check on the pdg code of the MC particle. if abs=kTRUE then check on the | |
3ff4a092 | 669 | //absolute value. |
563113d0 | 670 | // |
671 | TParticle* part = mcPart->Particle(); | |
672 | Int_t pdgCode = part->GetPdgCode(); | |
b95f6a36 | 673 | |
674 | if (abs) { | |
675 | pdgCode = TMath::Abs(pdgCode); | |
676 | pdg = TMath::Abs(pdg); | |
677 | } | |
678 | if (pdgCode != pdg ) return kFALSE; | |
679 | return kTRUE; | |
680 | } | |
681 | //______________________________ | |
682 | Bool_t AliCFParticleGenCuts::IsA(AliAODMCParticle *mcPart, Int_t pdg, Bool_t abs) { | |
683 | // | |
684 | //Check on the pdg code of the MC particle. if abs=kTRUE then check on the | |
3ff4a092 | 685 | //absolute value. |
b95f6a36 | 686 | // |
687 | Int_t pdgCode = mcPart->GetPdgCode(); | |
688 | ||
fc01457a | 689 | if (abs) { |
690 | pdgCode = TMath::Abs(pdgCode); | |
691 | pdg = TMath::Abs(pdg); | |
692 | } | |
107a3100 | 693 | if (pdgCode != pdg ) return kFALSE; |
563113d0 | 694 | return kTRUE; |
695 | } | |
696 | //______________________________ | |
107a3100 | 697 | void AliCFParticleGenCuts::SetEvtInfo(TObject* mcEvent) { |
563113d0 | 698 | // |
107a3100 | 699 | // Sets pointer to MC event information (AliMCEvent) |
563113d0 | 700 | // |
701 | ||
107a3100 | 702 | if (!mcEvent) { |
703 | AliError("Pointer to MC Event is null !"); | |
563113d0 | 704 | return; |
705 | } | |
706 | ||
107a3100 | 707 | TString className(mcEvent->ClassName()); |
b95f6a36 | 708 | if (className.CompareTo("AliMCEvent") != 0 && className.CompareTo("AliAODEvent") != 0) { |
709 | AliError("argument must point to an AliMCEvent or an AliAODEvent !"); | |
563113d0 | 710 | return ; |
711 | } | |
712 | ||
b95f6a36 | 713 | if (fIsAODMC) fMCInfo = dynamic_cast<AliAODEvent*>(mcEvent) ; |
714 | else fMCInfo = dynamic_cast<AliMCEvent*> (mcEvent) ; | |
563113d0 | 715 | } |