coverity defects
[u/mrichter/AliRoot.git] / CORRFW / AliCFParticleGenCuts.cxx
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
41 ClassImp(AliCFParticleGenCuts)
42
43 //______________________________
44 AliCFParticleGenCuts::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 //______________________________
90 AliCFParticleGenCuts::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 //______________________________
136 AliCFParticleGenCuts::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 //______________________________
182 AliCFParticleGenCuts& 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 //______________________________
233 Bool_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 //__________________________________________________________________________________
254 void 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 //__________________________________________________________________________________
426 void 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->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
538       }
539       if (!goodDecay) {
540         //second try inverting the order of the daughters
541         goodDecay = kTRUE ;
542         for (Int_t iDaughter = 0 ; iDaughter<nDaughters; iDaughter++) {
543           AliAODMCParticle* daug = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(1)-iDaughter));
544           if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
545         }
546       }
547     }
548     fCutValues->SetAt((Double32_t)goodDecay,kCutDecayChannel) ;
549   }
550   else fCutValues->SetAt((Double32_t)kTRUE,kCutDecayChannel);
551   
552   
553   // now array of cut is build, fill the bitmap consequently
554   Int_t iCutBit = -1;
555   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
556   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
557   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
558
559   ++iCutBit;
560   if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxXMin)
561     || ( fProdVtxRange2D && prodVtxXYmin >= 1))     fBitmap->SetBitNumber(iCutBit,kTRUE);
562
563   ++iCutBit;
564   if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxXMax)
565     || ( fProdVtxRange2D && prodVtxXYmax <= 1))     fBitmap->SetBitNumber(iCutBit,kTRUE);
566
567   ++iCutBit;
568   if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxYMin)
569     || ( fProdVtxRange2D && prodVtxXYmin >= 1))     fBitmap->SetBitNumber(iCutBit,kTRUE);
570
571   ++iCutBit;
572   if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxYMax)
573     || ( fProdVtxRange2D && prodVtxXYmax <= 1))     fBitmap->SetBitNumber(iCutBit,kTRUE);
574
575   if ( fCutValues->At(++iCutBit) > fProdVtxZMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
576   if ( fCutValues->At(++iCutBit) < fProdVtxZMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
577   if ( fCutValues->At(++iCutBit) > fDecayVtxXMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
578   if ( fCutValues->At(++iCutBit) < fDecayVtxXMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
579   if ( fCutValues->At(++iCutBit) > fDecayVtxYMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
580   if ( fCutValues->At(++iCutBit) < fDecayVtxYMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
581   if ( fCutValues->At(++iCutBit) > fDecayVtxZMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
582   if ( fCutValues->At(++iCutBit) < fDecayVtxZMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
583   if ( fCutValues->At(++iCutBit) > fDecayLengthMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
584   if ( fCutValues->At(++iCutBit) < fDecayLengthMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
585   if ( fCutValues->At(++iCutBit) > fDecayRxyMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
586   if ( fCutValues->At(++iCutBit) < fDecayRxyMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
587   if ( fCutValues->At(++iCutBit) != 0 )             fBitmap->SetBitNumber(iCutBit,kTRUE);
588 }
589
590
591 //__________________________________________________________________________________
592 void AliCFParticleGenCuts::FillHistograms(TObject* /*obj*/, Bool_t afterCuts)
593 {
594   //
595   // fill the QA histograms
596   //
597
598   for (int iCutNumber = 0; iCutNumber < kNCuts; iCutNumber++) 
599     fhQA[iCutNumber][afterCuts]->Fill(fCutValues->At(iCutNumber));
600
601   fhProdVtxXY[afterCuts]->Fill(fCutValues->At(4),fCutValues->At(5));
602
603   // fill cut statistics and cut correlation histograms with information from the bitmap
604   if (afterCuts) return;
605
606   // Number of single cuts in this class
607   UInt_t ncuts = fBitmap->GetNbits();
608   for(UInt_t bit=0; bit<ncuts;bit++) {
609     if (!fBitmap->TestBitNumber(bit)) {
610       fhCutStatistics->Fill(bit+1);
611       for (UInt_t bit2=bit; bit2<ncuts;bit2++) {
612         if (!fBitmap->TestBitNumber(bit2)) 
613           fhCutCorrelation->Fill(bit+1,bit2+1);
614       }
615     }
616   }
617 }
618
619 //__________________________________________________________________________________
620 void AliCFParticleGenCuts::AddQAHistograms(TList *qaList) {
621   //
622   // saves the histograms in a TList
623   //
624
625   DefineHistograms();
626
627   qaList->Add(fhCutStatistics);
628   qaList->Add(fhCutCorrelation);
629
630   for (Int_t j=0; j<kNStepQA; j++) {
631     qaList->Add(fhProdVtxXY[j]);
632     for(Int_t i=0; i<kNCuts; i++)
633       qaList->Add(fhQA[i][j]);
634   }
635 }
636
637 //__________________________________________________________________________________
638 void AliCFParticleGenCuts::DefineHistograms() {
639   //
640   // histograms for cut variables, cut statistics and cut correlations
641   //
642   Int_t color = 2;
643
644   // book cut statistics and cut correlation histograms
645   fhCutStatistics = new TH1F(Form("%s_cut_statistics",GetName()),"",kNCuts,0.5,kNCuts+0.5);
646   fhCutStatistics->SetLineWidth(2);
647   int k = 1;
648   fhCutStatistics->GetXaxis()->SetBinLabel(k,"charge")     ; k++;
649   fhCutStatistics->GetXaxis()->SetBinLabel(k,"prim/sec")   ; k++;
650   fhCutStatistics->GetXaxis()->SetBinLabel(k,"PDG")        ; k++;
651   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxXMin")    ; k++;
652   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxXMax")    ; k++;
653   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxYMin")    ; k++;
654   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxYMax")    ; k++;
655   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxZMin")    ; k++;
656   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMax")    ; k++;
657   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecXMin")    ; k++;
658   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecXMax")    ; k++;
659   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecYMin")    ; k++;
660   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecYMax")    ; k++;
661   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMin")    ; k++;
662   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMax")    ; k++;
663   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecLgthMin") ; k++;
664   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecLgthMax") ; k++;
665   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMin")  ; k++;
666   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMax")  ; k++;
667   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecChannel") ; k++;
668
669
670   fhCutCorrelation = new TH2F(Form("%s_cut_correlation",GetName()),"",kNCuts,0.5,kNCuts+0.5,kNCuts,0.5,kNCuts+0.5);
671   fhCutCorrelation->SetLineWidth(2);
672   for (k=1; k<=kNCuts; k++) {
673     fhCutCorrelation->GetXaxis()->SetBinLabel(k,fhCutStatistics->GetXaxis()->GetBinLabel(k));
674     fhCutCorrelation->GetYaxis()->SetBinLabel(k,fhCutStatistics->GetXaxis()->GetBinLabel(k));
675   }
676
677   Char_t str[5];
678   for (int i=0; i<kNStepQA; i++) {
679     if (i==0) snprintf(str,5," ");
680     else snprintf(str,5,"_cut");
681     fhQA[kCutCharge]      [i] = new TH1F(Form("%s_charge%s"      ,GetName(),str),"",2,0,2);
682     fhQA[kCutPrimSec]     [i] = new TH1F(Form("%s_primSec%s"     ,GetName(),str),"",2,0,2);
683     fhQA[kCutPDGCode]     [i] = new TH1F(Form("%s_pdgCode%s"     ,GetName(),str),"",2,0,2);
684     fhQA[kCutProdVtxXMin] [i] = new TH1F(Form("%s_prodVtxXMin%s" ,GetName(),str),"",100,-10,10);
685     fhQA[kCutProdVtxXMax] [i] = new TH1F(Form("%s_prodVtxXMax%s" ,GetName(),str),"",100,-10,10);
686     fhQA[kCutProdVtxYMin] [i] = new TH1F(Form("%s_prodVtxYMin%s" ,GetName(),str),"",100,-10,10);
687     fhQA[kCutProdVtxYMax] [i] = new TH1F(Form("%s_prodVtxYMax%s" ,GetName(),str),"",100,-10,10);
688     fhQA[kCutProdVtxZMin] [i] = new TH1F(Form("%s_prodVtxZMin%s" ,GetName(),str),"",100,-10,10);
689     fhQA[kCutProdVtxZMax] [i] = new TH1F(Form("%s_prodVtxZMax%s" ,GetName(),str),"",100,-10,10);
690     fhQA[kCutDecVtxXMin]  [i] = new TH1F(Form("%s_decVtxXMin%s"  ,GetName(),str),"",100,0,10);
691     fhQA[kCutDecVtxXMax]  [i] = new TH1F(Form("%s_decVtxXMax%s"  ,GetName(),str),"",100,0,10);
692     fhQA[kCutDecVtxYMin]  [i] = new TH1F(Form("%s_decVtxYMin%s"  ,GetName(),str),"",100,0,10);
693     fhQA[kCutDecVtxYMax]  [i] = new TH1F(Form("%s_decVtxYMax%s"  ,GetName(),str),"",100,0,10);
694     fhQA[kCutDecVtxZMin]  [i] = new TH1F(Form("%s_decVtxZMin%s"  ,GetName(),str),"",100,0,10);
695     fhQA[kCutDecVtxZMax]  [i] = new TH1F(Form("%s_decVtxZMax%s"  ,GetName(),str),"",100,0,10);
696     fhQA[kCutDecLgthMin]  [i] = new TH1F(Form("%s_decLengthMin%s",GetName(),str),"",100,0,10);
697     fhQA[kCutDecLgthMax]  [i] = new TH1F(Form("%s_decLengthMax%s",GetName(),str),"",100,0,10);
698     fhQA[kCutDecRxyMin]   [i] = new TH1F(Form("%s_decRxyMin%s"   ,GetName(),str),"",100,0,10);
699     fhQA[kCutDecRxyMax]   [i] = new TH1F(Form("%s_decRxyMax%s"   ,GetName(),str),"",100,0,10);
700     fhQA[kCutDecayChannel][i] = new TH1F(Form("%s_decayChannel%s",GetName(),str),"",2,0,2);
701     fhProdVtxXY           [i] = new TH2F(Form("%s_prodVtxXY%s"   ,GetName(),str),"",100,0,10,100,0,10);
702     fhProdVtxXY           [i] ->GetXaxis()->SetTitle("x_{production vertex}");
703     fhProdVtxXY           [i] ->GetYaxis()->SetTitle("y_{production vertex}");
704     fhQA[kCutProdVtxXMax] [i] ->GetXaxis()->SetTitle("x_{production vertex}");
705     fhQA[kCutProdVtxYMax] [i] ->GetXaxis()->SetTitle("y_{production vertex}");
706   }
707   for(Int_t i=0; i<kNCuts; i++) fhQA[i][1]->SetLineColor(color);
708 }
709
710
711 //______________________________
712 Bool_t AliCFParticleGenCuts::IsCharged(AliVParticle *mcPart) {
713   //
714   //check if particle is charged.
715   //
716   if (mcPart->Charge()==0) return kFALSE;
717   return kTRUE;
718 }
719 //______________________________
720 Bool_t AliCFParticleGenCuts::IsPrimary(AliMCParticle *mcPart) {
721   //
722   //check if particle is primary (standard definition)
723   //
724   
725   AliStack* stack = ((AliMCEvent*)fMCInfo)->Stack();
726
727   if (!stack->IsPhysicalPrimary(mcPart->GetLabel())) return kFALSE;
728   return kTRUE;
729 }
730 //______________________________
731 Bool_t AliCFParticleGenCuts::IsPrimary(AliAODMCParticle *mcPart) {
732   //
733   //check if particle is primary (standard definition)
734   //
735   
736   if (!mcPart->IsPhysicalPrimary()) return kFALSE;
737   return kTRUE;
738 }
739 //______________________________
740 Bool_t AliCFParticleGenCuts::IsPrimaryCharged(AliVParticle *mcPart) {
741   //
742   //check if a charged particle is primary (standard definition)
743   //
744
745   if (!fIsAODMC) {
746     if (!IsPrimary((AliMCParticle*)mcPart) || !IsCharged(mcPart)) return kFALSE ;
747   }
748   else {
749     if (!IsPrimary((AliAODMCParticle*)mcPart) || !IsCharged(mcPart)) return kFALSE ;
750   }
751   return kTRUE;
752 }
753 //______________________________
754 Bool_t AliCFParticleGenCuts::IsA(AliMCParticle *mcPart, Int_t pdg, Bool_t abs) {
755   //
756   //Check on the pdg code of the MC particle. if abs=kTRUE then check on the 
757   //absolute value. 
758   //
759   TParticle* part = mcPart->Particle();
760   Int_t pdgCode = part->GetPdgCode();
761
762   if (abs) {
763     pdgCode = TMath::Abs(pdgCode);
764     pdg = TMath::Abs(pdg);
765   }
766   if (pdgCode != pdg ) return kFALSE;
767   return kTRUE;
768 }
769 //______________________________
770 Bool_t AliCFParticleGenCuts::IsA(AliAODMCParticle *mcPart, Int_t pdg, Bool_t abs) {
771   //
772   //Check on the pdg code of the MC particle. if abs=kTRUE then check on the 
773   //absolute value. 
774   //
775   Int_t pdgCode = mcPart->GetPdgCode();
776   
777   if (abs) {
778     pdgCode = TMath::Abs(pdgCode);
779     pdg = TMath::Abs(pdg);
780   }
781   if (pdgCode != pdg ) return kFALSE;
782   return kTRUE;
783 }
784 //______________________________
785 void AliCFParticleGenCuts::SetMCEventInfo(const TObject* mcEvent) {
786   //
787   // Sets pointer to MC event information (AliMCEvent)
788   //
789
790   if (!mcEvent) {
791     AliError("Pointer to MC Event is null !");
792     return;
793   }
794   
795   TString className(mcEvent->ClassName());
796   if (className.CompareTo("AliMCEvent") != 0 && className.CompareTo("AliAODEvent") != 0) {
797     AliError("argument must point to an AliMCEvent or an AliAODEvent !");
798     return ;
799   }
800
801   fMCInfo = (AliVEvent*)mcEvent ;
802 }