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