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