13f25c8165446d468dbb56675658d69245e86091
[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     //now number of daughters is OK
318     if (goodDecay) {
319       // now check if decay channel is respected
320       // first try
321       for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) {
322         TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(iDaughter)) ;
323         if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
324       }
325       if (!goodDecay) {
326         //second try inverting the order of the daughters
327         goodDecay = kTRUE ;
328         for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) {
329           TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(nDaughters-(iDaughter+1))) ;
330           if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
331         }
332       }
333       if (!goodDecay && fRequireAbsolutePdg) {
334         //now tries inverting the sign of the daughters in case the anti-particle is also looked at
335         // third try
336         goodDecay = kTRUE ;
337         for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) {
338           TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(iDaughter)) ;
339           if (daug->GetPdgCode() != -fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
340         }
341         if (!goodDecay) {
342           //fourth try inverting the order of the daughters
343           goodDecay = kTRUE ;
344           for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) {
345             TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(nDaughters-(iDaughter+1))) ;
346             if (daug->GetPdgCode() != -fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
347           }
348         }
349       } //end check anti-particle
350     } //end # daughters OK
351     fCutValues->SetAt((Double32_t)goodDecay,kCutDecayChannel) ;
352   } //end require decay channel
353   else fCutValues->SetAt((Double32_t)kTRUE,kCutDecayChannel);
354   
355   
356   // now array of cut is build, fill the bitmap consequently
357   Int_t iCutBit = -1;
358   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
359   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
360   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
361   if ( fCutValues->At(++iCutBit) > fProdVtxXMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
362   if ( fCutValues->At(++iCutBit) < fProdVtxXMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
363   if ( fCutValues->At(++iCutBit) > fProdVtxYMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
364   if ( fCutValues->At(++iCutBit) < fProdVtxYMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
365   if ( fCutValues->At(++iCutBit) > fProdVtxZMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
366   if ( fCutValues->At(++iCutBit) < fProdVtxZMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
367   if ( fCutValues->At(++iCutBit) > fDecayVtxXMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
368   if ( fCutValues->At(++iCutBit) < fDecayVtxXMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
369   if ( fCutValues->At(++iCutBit) > fDecayVtxYMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
370   if ( fCutValues->At(++iCutBit) < fDecayVtxYMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
371   if ( fCutValues->At(++iCutBit) > fDecayVtxZMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
372   if ( fCutValues->At(++iCutBit) < fDecayVtxZMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
373   if ( fCutValues->At(++iCutBit) > fDecayLengthMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
374   if ( fCutValues->At(++iCutBit) < fDecayLengthMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
375   if ( fCutValues->At(++iCutBit) > fDecayRxyMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
376   if ( fCutValues->At(++iCutBit) < fDecayRxyMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
377   if ( fCutValues->At(++iCutBit) != 0 )             fBitmap->SetBitNumber(iCutBit,kTRUE);
378 }
379
380 //__________________________________________________________________________________
381 void AliCFParticleGenCuts::SelectionBitMap(AliAODMCParticle* mcPart)
382 {
383   //
384   // test if the track passes the single cuts
385   // and store the information in a bitmap
386   //
387   
388   for (UInt_t i=0; i<kNCuts; i++) {
389     fBitmap->SetBitNumber(i,kFALSE);
390     fCutValues->SetAt((Double32_t)0,i) ;
391   }
392
393   // fill the cut array
394   Double32_t partVx=(Double32_t)mcPart->Xv();
395   Double32_t partVy=(Double32_t)mcPart->Yv();
396   Double32_t partVz=(Double32_t)mcPart->Zv();
397
398   Double32_t decayVx=0.;
399   Double32_t decayVy=0.;
400   Double32_t decayVz=0.;
401   Double32_t decayL=0.;
402   Double32_t decayRxy=0.;
403
404   AliAODEvent* aod = dynamic_cast<AliAODEvent*>(fMCInfo);
405   TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
406   AliAODMCParticle* daughter = 0x0 ;
407
408   if (mcPart->GetDaughter(0)>0) {
409     daughter = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(0)));
410
411     decayVx=(Double32_t)daughter->Xv();
412     decayVy=(Double32_t)daughter->Yv();
413     decayVz=(Double32_t)daughter->Zv();
414     decayL = TMath::Sqrt(TMath::Power(partVx-decayVx,2) + 
415                          TMath::Power(partVy-decayVy,2) + 
416                          TMath::Power(partVz-decayVz,2) ) ;
417     decayRxy = TMath::Sqrt(TMath::Power(decayVx,2) + TMath::Power(decayVy,2) ) ;
418     
419   }
420   
421   fCutValues->SetAt(partVx  ,kCutProdVtxXMin) ;
422   fCutValues->SetAt(partVy  ,kCutProdVtxYMin) ;
423   fCutValues->SetAt(partVz  ,kCutProdVtxZMin) ;
424   fCutValues->SetAt(partVx  ,kCutProdVtxXMax) ;
425   fCutValues->SetAt(partVy  ,kCutProdVtxYMax) ;
426   fCutValues->SetAt(partVz  ,kCutProdVtxZMax) ;
427   fCutValues->SetAt(decayVx ,kCutDecVtxXMin)  ;
428   fCutValues->SetAt(decayVy ,kCutDecVtxYMin)  ;
429   fCutValues->SetAt(decayVz ,kCutDecVtxZMin)  ;
430   fCutValues->SetAt(decayVx ,kCutDecVtxXMax)  ;
431   fCutValues->SetAt(decayVy ,kCutDecVtxYMax)  ;
432   fCutValues->SetAt(decayVz ,kCutDecVtxZMax)  ;
433   fCutValues->SetAt(decayL  ,kCutDecLgthMin)  ;
434   fCutValues->SetAt(decayL  ,kCutDecLgthMax)  ;
435   fCutValues->SetAt(decayRxy,kCutDecRxyMin)   ;
436   fCutValues->SetAt(decayRxy,kCutDecRxyMax)   ;
437   
438   // cut on charge
439   if ( fRequireIsCharged || fRequireIsNeutral ) {
440     if (fRequireIsCharged &&  IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
441     if (fRequireIsNeutral && !IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
442   } 
443   else fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
444   
445   // cut on primary/secondary
446   if ( fRequireIsPrimary || fRequireIsSecondary) {
447     if (fRequireIsPrimary   &&  IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
448     if (fRequireIsSecondary && !IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
449   } 
450   else fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
451   
452   // cut on PDG code
453   if ( fRequirePdgCode ) {
454     if (IsA(mcPart,fPdgCode,fRequireAbsolutePdg)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
455   }
456   else fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
457   
458   // cut on decay channel
459   if ( fDecayChannel ) {
460     Bool_t goodDecay = kTRUE ;
461     Short_t nDaughters = 0 ;
462     if (mcPart->GetDaughter(0) >=0) nDaughters = 1 + mcPart->GetDaughter(1) - mcPart->GetDaughter(0);
463     
464     if (nDaughters != fDecayChannel->NDaughters()) goodDecay = kFALSE ;
465     if (goodDecay) {
466       // now check if decay channel is respected
467       // first try
468       for (Int_t iDaughter = 0 ; iDaughter<nDaughters; iDaughter++) {
469         AliAODMCParticle* daug = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(0)+iDaughter));
470         if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
471       }
472       if (!goodDecay) {
473         //second try inverting the order of the daughters
474         goodDecay = kTRUE ;
475         for (Int_t iDaughter = 0 ; iDaughter<nDaughters; iDaughter++) {
476           AliAODMCParticle* daug = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(1)-iDaughter));
477           if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
478         }
479       }
480     }
481     fCutValues->SetAt((Double32_t)goodDecay,kCutDecayChannel) ;
482   }
483   else fCutValues->SetAt((Double32_t)kTRUE,kCutDecayChannel);
484   
485   
486   // now array of cut is build, fill the bitmap consequently
487   Int_t iCutBit = -1;
488   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
489   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
490   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
491   if ( fCutValues->At(++iCutBit) > fProdVtxXMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
492   if ( fCutValues->At(++iCutBit) < fProdVtxXMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
493   if ( fCutValues->At(++iCutBit) > fProdVtxYMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
494   if ( fCutValues->At(++iCutBit) < fProdVtxYMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
495   if ( fCutValues->At(++iCutBit) > fProdVtxZMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
496   if ( fCutValues->At(++iCutBit) < fProdVtxZMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
497   if ( fCutValues->At(++iCutBit) > fDecayVtxXMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
498   if ( fCutValues->At(++iCutBit) < fDecayVtxXMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
499   if ( fCutValues->At(++iCutBit) > fDecayVtxYMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
500   if ( fCutValues->At(++iCutBit) < fDecayVtxYMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
501   if ( fCutValues->At(++iCutBit) > fDecayVtxZMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
502   if ( fCutValues->At(++iCutBit) < fDecayVtxZMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
503   if ( fCutValues->At(++iCutBit) > fDecayLengthMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
504   if ( fCutValues->At(++iCutBit) < fDecayLengthMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
505   if ( fCutValues->At(++iCutBit) > fDecayRxyMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
506   if ( fCutValues->At(++iCutBit) < fDecayRxyMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
507   if ( fCutValues->At(++iCutBit) != 0 )             fBitmap->SetBitNumber(iCutBit,kTRUE);
508 }
509
510
511 //__________________________________________________________________________________
512 void AliCFParticleGenCuts::FillHistograms(TObject* /*obj*/, Bool_t afterCuts)
513 {
514   //
515   // fill the QA histograms
516   //
517
518   for (int iCutNumber = 0; iCutNumber < kNCuts; iCutNumber++) 
519     fhQA[iCutNumber][afterCuts]->Fill(fCutValues->At(iCutNumber));
520
521   // fill cut statistics and cut correlation histograms with information from the bitmap
522   if (afterCuts) return;
523
524   // Number of single cuts in this class
525   UInt_t ncuts = fBitmap->GetNbits();
526   for(UInt_t bit=0; bit<ncuts;bit++) {
527     if (!fBitmap->TestBitNumber(bit)) {
528       fhCutStatistics->Fill(bit+1);
529       for (UInt_t bit2=bit; bit2<ncuts;bit2++) {
530         if (!fBitmap->TestBitNumber(bit2)) 
531           fhCutCorrelation->Fill(bit+1,bit2+1);
532       }
533     }
534   }
535 }
536
537 //__________________________________________________________________________________
538 void AliCFParticleGenCuts::AddQAHistograms(TList *qaList) {
539   //
540   // saves the histograms in a TList
541   //
542
543   DefineHistograms();
544
545   qaList->Add(fhCutStatistics);
546   qaList->Add(fhCutCorrelation);
547
548   for (Int_t j=0; j<kNStepQA; j++) {
549     for(Int_t i=0; i<kNCuts; i++)
550       qaList->Add(fhQA[i][j]);
551   }
552 }
553
554 //__________________________________________________________________________________
555 void AliCFParticleGenCuts::DefineHistograms() {
556   //
557   // histograms for cut variables, cut statistics and cut correlations
558   //
559   Int_t color = 2;
560
561   // book cut statistics and cut correlation histograms
562   fhCutStatistics = new TH1F(Form("%s_cut_statistics",GetName()),"",kNCuts,0.5,kNCuts+0.5);
563   fhCutStatistics->SetLineWidth(2);
564   int k = 1;
565   fhCutStatistics->GetXaxis()->SetBinLabel(k,"charge")     ; k++;
566   fhCutStatistics->GetXaxis()->SetBinLabel(k,"prim/sec")   ; k++;
567   fhCutStatistics->GetXaxis()->SetBinLabel(k,"PDG")        ; k++;
568   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxXMin")    ; k++;
569   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxXMax")    ; k++;
570   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxYMin")    ; k++;
571   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxYMax")    ; k++;
572   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxZMin")    ; k++;
573   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMax")    ; k++;
574   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecXMin")    ; k++;
575   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecXMax")    ; k++;
576   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecYMin")    ; k++;
577   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecYMax")    ; k++;
578   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMin")    ; k++;
579   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMax")    ; k++;
580   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecLgthMin") ; k++;
581   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecLgthMax") ; k++;
582   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMin")  ; k++;
583   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMax")  ; k++;
584   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecChannel") ; k++;
585
586
587   fhCutCorrelation = new TH2F(Form("%s_cut_correlation",GetName()),"",kNCuts,0.5,kNCuts+0.5,kNCuts,0.5,kNCuts+0.5);
588   fhCutCorrelation->SetLineWidth(2);
589   for (k=1; k<=kNCuts; k++) {
590     fhCutCorrelation->GetXaxis()->SetBinLabel(k,fhCutStatistics->GetXaxis()->GetBinLabel(k));
591     fhCutCorrelation->GetYaxis()->SetBinLabel(k,fhCutStatistics->GetXaxis()->GetBinLabel(k));
592   }
593
594   Char_t str[256];
595   for (int i=0; i<kNStepQA; i++) {
596     if (i==0) sprintf(str," ");
597     else sprintf(str,"_cut");
598     fhQA[kCutCharge]      [i] = new TH1F(Form("%s_charge%s"      ,GetName(),str),"",2,0,2);
599     fhQA[kCutPrimSec]     [i] = new TH1F(Form("%s_primSec%s"     ,GetName(),str),"",2,0,2);
600     fhQA[kCutPDGCode]     [i] = new TH1F(Form("%s_pdgCode%s"     ,GetName(),str),"",2,0,2);
601     fhQA[kCutProdVtxXMin] [i] = new TH1F(Form("%s_prodVtxXMin%s" ,GetName(),str),"",100,0,10);
602     fhQA[kCutProdVtxXMax] [i] = new TH1F(Form("%s_prodVtxXMax%s" ,GetName(),str),"",100,0,10);
603     fhQA[kCutProdVtxYMin] [i] = new TH1F(Form("%s_prodVtxYMin%s" ,GetName(),str),"",100,0,10);
604     fhQA[kCutProdVtxYMax] [i] = new TH1F(Form("%s_prodVtxYMax%s" ,GetName(),str),"",100,0,10);
605     fhQA[kCutProdVtxZMin] [i] = new TH1F(Form("%s_prodVtxZMin%s" ,GetName(),str),"",100,0,10);
606     fhQA[kCutProdVtxZMax] [i] = new TH1F(Form("%s_prodVtxZMax%s" ,GetName(),str),"",100,0,10);
607     fhQA[kCutDecVtxXMin]  [i] = new TH1F(Form("%s_decVtxXMin%s"  ,GetName(),str),"",100,0,10);
608     fhQA[kCutDecVtxXMax]  [i] = new TH1F(Form("%s_decVtxXMax%s"  ,GetName(),str),"",100,0,10);
609     fhQA[kCutDecVtxYMin]  [i] = new TH1F(Form("%s_decVtxYMin%s"  ,GetName(),str),"",100,0,10);
610     fhQA[kCutDecVtxYMax]  [i] = new TH1F(Form("%s_decVtxYMax%s"  ,GetName(),str),"",100,0,10);
611     fhQA[kCutDecVtxZMin]  [i] = new TH1F(Form("%s_decVtxZMin%s"  ,GetName(),str),"",100,0,10);
612     fhQA[kCutDecVtxZMax]  [i] = new TH1F(Form("%s_decVtxZMax%s"  ,GetName(),str),"",100,0,10);
613     fhQA[kCutDecLgthMin]  [i] = new TH1F(Form("%s_decLengthMin%s",GetName(),str),"",100,0,10);
614     fhQA[kCutDecLgthMax]  [i] = new TH1F(Form("%s_decLengthMax%s",GetName(),str),"",100,0,10);
615     fhQA[kCutDecRxyMin]   [i] = new TH1F(Form("%s_decRxyMin%s"   ,GetName(),str),"",100,0,10);
616     fhQA[kCutDecRxyMax]   [i] = new TH1F(Form("%s_decRxyMax%s"   ,GetName(),str),"",100,0,10);
617     fhQA[kCutDecayChannel][i] = new TH1F(Form("%s_decayChannel%s",GetName(),str),"",2,0,2);
618   }
619   for(Int_t i=0; i<kNCuts; i++) fhQA[i][1]->SetLineColor(color);
620 }
621
622
623 //______________________________
624 Bool_t AliCFParticleGenCuts::IsCharged(AliVParticle *mcPart) {
625   //
626   //check if particle is charged.
627   //
628   if (mcPart->Charge()==0) return kFALSE;
629   return kTRUE;
630 }
631 //______________________________
632 Bool_t AliCFParticleGenCuts::IsPrimary(AliMCParticle *mcPart) {
633   //
634   //check if particle is primary (standard definition)
635   //
636   
637   AliStack* stack = ((AliMCEvent*)fMCInfo)->Stack();
638
639   if (!stack->IsPhysicalPrimary(mcPart->GetLabel())) return kFALSE;
640   return kTRUE;
641 }
642 //______________________________
643 Bool_t AliCFParticleGenCuts::IsPrimary(AliAODMCParticle *mcPart) {
644   //
645   //check if particle is primary (standard definition)
646   //
647   
648   if (!mcPart->IsPhysicalPrimary()) return kFALSE;
649   return kTRUE;
650 }
651 //______________________________
652 Bool_t AliCFParticleGenCuts::IsPrimaryCharged(AliVParticle *mcPart) {
653   //
654   //check if a charged particle is primary (standard definition)
655   //
656
657   if (!fIsAODMC) {
658     if (!IsPrimary((AliMCParticle*)mcPart) || !IsCharged(mcPart)) return kFALSE ;
659   }
660   else {
661     if (!IsPrimary((AliAODMCParticle*)mcPart) || !IsCharged(mcPart)) return kFALSE ;
662   }
663   return kTRUE;
664 }
665 //______________________________
666 Bool_t AliCFParticleGenCuts::IsA(AliMCParticle *mcPart, Int_t pdg, Bool_t abs) {
667   //
668   //Check on the pdg code of the MC particle. if abs=kTRUE then check on the 
669   //absolute value. 
670   //
671   TParticle* part = mcPart->Particle();
672   Int_t pdgCode = part->GetPdgCode();
673
674   if (abs) {
675     pdgCode = TMath::Abs(pdgCode);
676     pdg = TMath::Abs(pdg);
677   }
678   if (pdgCode != pdg ) return kFALSE;
679   return kTRUE;
680 }
681 //______________________________
682 Bool_t AliCFParticleGenCuts::IsA(AliAODMCParticle *mcPart, Int_t pdg, Bool_t abs) {
683   //
684   //Check on the pdg code of the MC particle. if abs=kTRUE then check on the 
685   //absolute value. 
686   //
687   Int_t pdgCode = mcPart->GetPdgCode();
688   
689   if (abs) {
690     pdgCode = TMath::Abs(pdgCode);
691     pdg = TMath::Abs(pdg);
692   }
693   if (pdgCode != pdg ) return kFALSE;
694   return kTRUE;
695 }
696 //______________________________
697 void AliCFParticleGenCuts::SetEvtInfo(TObject* mcEvent) {
698   //
699   // Sets pointer to MC event information (AliMCEvent)
700   //
701
702   if (!mcEvent) {
703     AliError("Pointer to MC Event is null !");
704     return;
705   }
706   
707   TString className(mcEvent->ClassName());
708   if (className.CompareTo("AliMCEvent") != 0 && className.CompareTo("AliAODEvent") != 0) {
709     AliError("argument must point to an AliMCEvent or an AliAODEvent !");
710     return ;
711   }
712   
713   if (fIsAODMC) fMCInfo = dynamic_cast<AliAODEvent*>(mcEvent) ;
714   else          fMCInfo = dynamic_cast<AliMCEvent*> (mcEvent) ;
715 }