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