Removing duplicate files
[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
39 ClassImp(AliCFParticleGenCuts)
40
41 //______________________________
42 AliCFParticleGenCuts::AliCFParticleGenCuts() : 
43   AliCFCutBase(),
44   fMCInfo(0x0),
45   fRequireIsCharged(0),
46   fRequireIsNeutral(0),
47   fRequireIsPrimary(0),
48   fRequireIsSecondary(0),
49   fRequirePdgCode(0),
50   fPdgCode(0),
51   fProdVtxXMin (-1.e+09),
52   fProdVtxYMin (-1.e+09),
53   fProdVtxZMin (-1.e+09),
54   fProdVtxXMax ( 1.e+09),
55   fProdVtxYMax ( 1.e+09),
56   fProdVtxZMax ( 1.e+09),
57   fDecayVtxXMin(-1.e+09),
58   fDecayVtxYMin(-1.e+09),
59   fDecayVtxZMin(-1.e+09),
60   fDecayVtxXMax( 1.e+09),
61   fDecayVtxYMax( 1.e+09),
62   fDecayVtxZMax( 1.e+09),
63   fDecayLengthMin(-1.),
64   fDecayLengthMax(1.e+09),
65   fDecayRxyMin(-1),
66   fDecayRxyMax(1.e+09),
67   fDecayChannel(0x0),
68   fhCutStatistics(0x0),
69   fhCutCorrelation(0x0),
70   fCutValues(new TArrayF(kNCuts)),
71   fBitmap(new TBits(0))
72 {
73   //
74   //ctor
75   //
76   for (int i=0; i<kNCuts; i++) 
77     for (int j=0; j<kNStepQA; j++) 
78       fhQA[i][j]=0x0;
79 }
80
81 //______________________________
82 AliCFParticleGenCuts::AliCFParticleGenCuts(const Char_t* name, const Char_t* title) : 
83   AliCFCutBase(name,title),
84   fMCInfo(0x0),
85   fRequireIsCharged(0),
86   fRequireIsNeutral(0),
87   fRequireIsPrimary(0),
88   fRequireIsSecondary(0),
89   fRequirePdgCode(0),
90   fPdgCode(0),
91   fProdVtxXMin (-1.e+09),
92   fProdVtxYMin (-1.e+09),
93   fProdVtxZMin (-1.e+09),
94   fProdVtxXMax ( 1.e+09),
95   fProdVtxYMax ( 1.e+09),
96   fProdVtxZMax ( 1.e+09),
97   fDecayVtxXMin(-1.e+09),
98   fDecayVtxYMin(-1.e+09),
99   fDecayVtxZMin(-1.e+09),
100   fDecayVtxXMax( 1.e+09),
101   fDecayVtxYMax( 1.e+09),
102   fDecayVtxZMax( 1.e+09),
103   fDecayLengthMin(-1.),
104   fDecayLengthMax(1.e+09),
105   fDecayRxyMin(-1.),
106   fDecayRxyMax(1.e+09),
107   fDecayChannel(0x0),
108   fhCutStatistics(0x0),
109   fhCutCorrelation(0x0),
110   fCutValues(new TArrayF(kNCuts)),
111   fBitmap(new TBits(0))
112 {
113   //
114   //ctor
115   //
116   for (int i=0; i<kNCuts; i++) 
117     for (int j=0; j<kNStepQA; j++) 
118       fhQA[i][j]=0x0;
119 }
120
121 //______________________________
122 AliCFParticleGenCuts::AliCFParticleGenCuts(const AliCFParticleGenCuts& c) : 
123   AliCFCutBase(c),
124   fMCInfo(c.fMCInfo),
125   fRequireIsCharged(c.fRequireIsCharged),
126   fRequireIsNeutral(c.fRequireIsNeutral),
127   fRequireIsPrimary(c.fRequireIsPrimary),
128   fRequireIsSecondary(c.fRequireIsSecondary),
129   fRequirePdgCode(c.fRequirePdgCode),
130   fPdgCode(c.fPdgCode),
131   fProdVtxXMin (c.fProdVtxXMin),
132   fProdVtxYMin (c.fProdVtxYMin),
133   fProdVtxZMin (c.fProdVtxZMin),
134   fProdVtxXMax (c.fProdVtxXMax),
135   fProdVtxYMax (c.fProdVtxYMax),
136   fProdVtxZMax (c.fProdVtxZMax),
137   fDecayVtxXMin(c.fDecayVtxXMin),
138   fDecayVtxYMin(c.fDecayVtxYMin),
139   fDecayVtxZMin(c.fDecayVtxZMin),
140   fDecayVtxXMax(c.fDecayVtxXMax),
141   fDecayVtxYMax(c.fDecayVtxYMax),
142   fDecayVtxZMax(c.fDecayVtxZMax),
143   fDecayLengthMin(c.fDecayLengthMin),
144   fDecayLengthMax(c.fDecayLengthMin),
145   fDecayRxyMin(c.fDecayLengthMin),
146   fDecayRxyMax(c.fDecayLengthMin),
147   fDecayChannel(c.fDecayChannel),
148   fhCutStatistics(new TH1F(*c.fhCutStatistics)),
149   fhCutCorrelation(new TH2F(*c.fhCutCorrelation)),
150   fCutValues(new TArrayF(*c.fCutValues)),
151   fBitmap(new TBits(*c.fBitmap))
152 {
153   //
154   //copy ctor
155   //
156   for (int i=0; i<kNCuts; i++)
157     for (int j=0; j<kNStepQA; j++)
158       fhQA[i][j]=(TH1F*)c.fhQA[i][j]->Clone();
159 }
160
161 //______________________________
162 AliCFParticleGenCuts& AliCFParticleGenCuts::operator=(const AliCFParticleGenCuts& c)
163 {
164   //
165   // Assignment operator
166   //
167   if (this != &c) {
168     AliCFCutBase::operator=(c) ;
169     fMCInfo=c.fMCInfo;
170     fRequireIsCharged=c.fRequireIsCharged;
171     fRequireIsNeutral=c.fRequireIsNeutral;
172     fRequireIsPrimary=c.fRequireIsPrimary;
173     fRequireIsSecondary=c.fRequireIsSecondary;
174     fRequirePdgCode=c.fRequirePdgCode;
175     fPdgCode=c.fPdgCode;
176     fProdVtxXMin=c.fProdVtxXMin;
177     fProdVtxYMin=c.fProdVtxYMin;
178     fProdVtxZMin=c.fProdVtxZMin;
179     fProdVtxXMax=c.fProdVtxXMax;
180     fProdVtxYMax=c.fProdVtxYMax;
181     fProdVtxZMax=c.fProdVtxZMax;
182     fDecayVtxXMin=c.fDecayVtxXMin;
183     fDecayVtxYMin=c.fDecayVtxYMin;
184     fDecayVtxZMin=c.fDecayVtxZMin;
185     fDecayVtxXMax=c.fDecayVtxXMax;
186     fDecayVtxYMax=c.fDecayVtxYMax;
187     fDecayVtxZMax=c.fDecayVtxZMax;      
188     fDecayLengthMin=c.fDecayVtxZMax;
189     fDecayLengthMax=c.fDecayLengthMax;
190     fDecayRxyMin=c.fDecayRxyMin;
191     fDecayRxyMax=c.fDecayRxyMax;
192     fDecayChannel=c.fDecayChannel;
193     fCutValues=new TArrayF(*c.fCutValues);
194     fBitmap=new TBits(*c.fBitmap);
195     
196     if (fhCutStatistics)  fhCutStatistics =new TH1F(*c.fhCutStatistics) ;
197     if (fhCutCorrelation) fhCutCorrelation=new TH2F(*c.fhCutCorrelation);
198     
199     for (int i=0; i<kNCuts; i++)
200       for (int j=0; j<kNStepQA; j++)
201         fhQA[i][j]=(TH1F*)c.fhQA[i][j]->Clone();
202   }
203   return *this ;
204 }
205
206 //______________________________
207 Bool_t AliCFParticleGenCuts::IsSelected(TObject* obj) {
208   //
209   // check if selections on 'obj' are passed
210   // 'obj' must be an AliMCParticle
211   //
212   
213   SelectionBitMap(obj);
214
215   if (fIsQAOn) FillHistograms(obj,0);
216
217   for (UInt_t icut=0; icut<fBitmap->GetNbits();icut++)
218     if (!fBitmap->TestBitNumber(icut)) return kFALSE ; 
219   
220   if (fIsQAOn) FillHistograms(obj,1);
221   return kTRUE;
222 }
223              
224 //__________________________________________________________________________________
225 void AliCFParticleGenCuts::SelectionBitMap(TObject* obj)
226 {
227   //
228   // test if the track passes the single cuts
229   // and store the information in a bitmap
230   //
231   
232   for (UInt_t i=0; i<kNCuts; i++) {
233     fBitmap->SetBitNumber(i,kFALSE);
234     fCutValues->SetAt((Double32_t)0,i) ;
235   }
236
237   if (!obj) return  ;
238   TString className(obj->ClassName());
239   if (className.CompareTo("AliMCParticle") != 0) {
240     AliError("argument must point to an AliMCParticle !");
241     return ;
242   }
243   
244   AliMCParticle* mcPart = dynamic_cast<AliMCParticle*>(obj) ;
245   TParticle* part = mcPart->Particle();
246   AliStack*  stack = fMCInfo->Stack();
247
248
249   // fill the cut array
250   Double32_t partVx=(Double32_t)part->Vx();
251   Double32_t partVy=(Double32_t)part->Vy();
252   Double32_t partVz=(Double32_t)part->Vz();
253
254   TParticle* daughter=0x0;
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   if ( part->GetNDaughters() > 0 ) {
262     daughter = stack->Particle(part->GetFirstDaughter()) ;
263     decayVx=(Double32_t)daughter->Vx();
264     decayVy=(Double32_t)daughter->Vy();
265     decayVz=(Double32_t)daughter->Vz();
266     decayL = TMath::Sqrt(TMath::Power(partVx-decayVx,2) + 
267                          TMath::Power(partVy-decayVy,2) + 
268                          TMath::Power(partVz-decayVz,2) ) ;
269     decayRxy = TMath::Sqrt(TMath::Power(decayVx,2) + TMath::Power(decayVy,2) ) ;
270   }
271
272   fCutValues->SetAt(partVx  ,kCutProdVtxXMin) ;
273   fCutValues->SetAt(partVy  ,kCutProdVtxYMin) ;
274   fCutValues->SetAt(partVz  ,kCutProdVtxZMin) ;
275   fCutValues->SetAt(partVx  ,kCutProdVtxXMax) ;
276   fCutValues->SetAt(partVy  ,kCutProdVtxYMax) ;
277   fCutValues->SetAt(partVz  ,kCutProdVtxZMax) ;
278   fCutValues->SetAt(decayVx ,kCutDecVtxXMin)  ;
279   fCutValues->SetAt(decayVy ,kCutDecVtxYMin)  ;
280   fCutValues->SetAt(decayVz ,kCutDecVtxZMin)  ;
281   fCutValues->SetAt(decayVx ,kCutDecVtxXMax)  ;
282   fCutValues->SetAt(decayVy ,kCutDecVtxYMax)  ;
283   fCutValues->SetAt(decayVz ,kCutDecVtxZMax)  ;
284   fCutValues->SetAt(decayL  ,kCutDecLgthMin)  ;
285   fCutValues->SetAt(decayL  ,kCutDecLgthMax)  ;
286   fCutValues->SetAt(decayRxy,kCutDecRxyMin)   ;
287   fCutValues->SetAt(decayRxy,kCutDecRxyMax)   ;
288   
289   // cut on charge
290   if ( fRequireIsCharged || fRequireIsNeutral ) {
291     if (fRequireIsCharged &&  IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
292     if (fRequireIsNeutral && !IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
293   } 
294   else fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
295   
296   // cut on primary/secondary
297   if ( fRequireIsPrimary || fRequireIsSecondary) {
298     if (fRequireIsPrimary   &&  IsPrimary(mcPart,stack)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
299     if (fRequireIsSecondary && !IsPrimary(mcPart,stack)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
300   } 
301   else fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
302   
303   // cut on PDG code
304   if ( fRequirePdgCode ) {
305     if (IsA(mcPart,fPdgCode,kFALSE)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
306   }
307   else fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
308   
309   // cut on decay channel
310   if ( fDecayChannel ) {
311     Bool_t goodDecay = kTRUE ;
312     Short_t nDaughters = mcPart->Particle()->GetNDaughters() ;
313     if (nDaughters != fDecayChannel->NDaughters()) goodDecay = kFALSE ;
314     if (goodDecay) {
315       // now check if decay channel is respected
316       // first try
317       for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) {
318         TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(iDaughter)) ;
319         if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
320       }
321       if (!goodDecay) {
322         //second try inverting the order of the daughters
323         goodDecay = kTRUE ;
324         for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) {
325           TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(nDaughters-(iDaughter+1))) ;
326           if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
327         }
328       }
329     }
330     fCutValues->SetAt((Double32_t)goodDecay,kCutDecayChannel) ;
331   }
332   else fCutValues->SetAt((Double32_t)kTRUE,kCutDecayChannel);
333   
334   
335   // now array of cut is build, fill the bitmap consequently
336   Int_t iCutBit = -1;
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) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
340   if ( fCutValues->At(++iCutBit) > fProdVtxXMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
341   if ( fCutValues->At(++iCutBit) < fProdVtxXMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
342   if ( fCutValues->At(++iCutBit) > fProdVtxYMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
343   if ( fCutValues->At(++iCutBit) < fProdVtxYMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
344   if ( fCutValues->At(++iCutBit) > fProdVtxZMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
345   if ( fCutValues->At(++iCutBit) < fProdVtxZMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
346   if ( fCutValues->At(++iCutBit) > fDecayVtxXMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
347   if ( fCutValues->At(++iCutBit) < fDecayVtxXMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
348   if ( fCutValues->At(++iCutBit) > fDecayVtxYMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
349   if ( fCutValues->At(++iCutBit) < fDecayVtxYMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
350   if ( fCutValues->At(++iCutBit) > fDecayVtxZMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
351   if ( fCutValues->At(++iCutBit) < fDecayVtxZMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
352   if ( fCutValues->At(++iCutBit) > fDecayLengthMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
353   if ( fCutValues->At(++iCutBit) < fDecayLengthMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
354   if ( fCutValues->At(++iCutBit) > fDecayRxyMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
355   if ( fCutValues->At(++iCutBit) < fDecayRxyMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
356   if ( fCutValues->At(++iCutBit) != 0 )             fBitmap->SetBitNumber(iCutBit,kTRUE);
357 }
358
359
360 //__________________________________________________________________________________
361 void AliCFParticleGenCuts::FillHistograms(TObject* /*obj*/, Bool_t afterCuts)
362 {
363   //
364   // fill the QA histograms
365   //
366
367   for (int iCutNumber = 0; iCutNumber < kNCuts; iCutNumber++) 
368     fhQA[iCutNumber][afterCuts]->Fill(fCutValues->At(iCutNumber));
369
370   // fill cut statistics and cut correlation histograms with information from the bitmap
371   if (afterCuts) return;
372
373   // Number of single cuts in this class
374   UInt_t ncuts = fBitmap->GetNbits();
375   for(UInt_t bit=0; bit<ncuts;bit++) {
376     if (!fBitmap->TestBitNumber(bit)) {
377       fhCutStatistics->Fill(bit+1);
378       for (UInt_t bit2=bit; bit2<ncuts;bit2++) {
379         if (!fBitmap->TestBitNumber(bit2)) 
380           fhCutCorrelation->Fill(bit+1,bit2+1);
381       }
382     }
383   }
384 }
385
386 //__________________________________________________________________________________
387 void AliCFParticleGenCuts::AddQAHistograms(TList *qaList) {
388   //
389   // saves the histograms in a TList
390   //
391
392   DefineHistograms();
393
394   qaList->Add(fhCutStatistics);
395   qaList->Add(fhCutCorrelation);
396
397   for (Int_t j=0; j<kNStepQA; j++) {
398     for(Int_t i=0; i<kNCuts; i++)
399       qaList->Add(fhQA[i][j]);
400   }
401 }
402
403 //__________________________________________________________________________________
404 void AliCFParticleGenCuts::DefineHistograms() {
405   //
406   // histograms for cut variables, cut statistics and cut correlations
407   //
408   Int_t color = 2;
409
410   // book cut statistics and cut correlation histograms
411   fhCutStatistics = new TH1F(Form("%s_cut_statistics",GetName()),"",kNCuts,0.5,kNCuts+0.5);
412   fhCutStatistics->SetLineWidth(2);
413   int k = 1;
414   fhCutStatistics->GetXaxis()->SetBinLabel(k,"charge")     ; k++;
415   fhCutStatistics->GetXaxis()->SetBinLabel(k,"prim/sec")   ; k++;
416   fhCutStatistics->GetXaxis()->SetBinLabel(k,"PDG")        ; k++;
417   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxXMin")    ; k++;
418   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxXMax")    ; k++;
419   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxYMin")    ; k++;
420   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxYMax")    ; k++;
421   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxZMin")    ; k++;
422   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMax")    ; k++;
423   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecXMin")    ; k++;
424   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecXMax")    ; k++;
425   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecYMin")    ; k++;
426   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecYMax")    ; k++;
427   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMin")    ; k++;
428   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMax")    ; k++;
429   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecLgthMin") ; k++;
430   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecLgthMax") ; k++;
431   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMin")  ; k++;
432   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMax")  ; k++;
433   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecChannel") ; k++;
434
435
436   fhCutCorrelation = new TH2F(Form("%s_cut_correlation",GetName()),"",kNCuts,0.5,kNCuts+0.5,kNCuts,0.5,kNCuts+0.5);
437   fhCutCorrelation->SetLineWidth(2);
438   for (k=1; k<=kNCuts; k++) {
439     fhCutCorrelation->GetXaxis()->SetBinLabel(k,fhCutStatistics->GetXaxis()->GetBinLabel(k));
440     fhCutCorrelation->GetYaxis()->SetBinLabel(k,fhCutStatistics->GetXaxis()->GetBinLabel(k));
441   }
442
443   Char_t str[256];
444   for (int i=0; i<kNStepQA; i++) {
445     if (i==0) sprintf(str," ");
446     else sprintf(str,"_cut");
447     fhQA[kCutCharge]      [i] = new TH1F(Form("%s_charge%s"      ,GetName(),str),"",2,0,2);
448     fhQA[kCutPrimSec]     [i] = new TH1F(Form("%s_primSec%s"     ,GetName(),str),"",2,0,2);
449     fhQA[kCutPDGCode]     [i] = new TH1F(Form("%s_pdgCode%s"     ,GetName(),str),"",2,0,2);
450     fhQA[kCutProdVtxXMin] [i] = new TH1F(Form("%s_prodVtxXMin%s" ,GetName(),str),"",100,0,10);
451     fhQA[kCutProdVtxXMax] [i] = new TH1F(Form("%s_prodVtxXMax%s" ,GetName(),str),"",100,0,10);
452     fhQA[kCutProdVtxYMin] [i] = new TH1F(Form("%s_prodVtxYMin%s" ,GetName(),str),"",100,0,10);
453     fhQA[kCutProdVtxYMax] [i] = new TH1F(Form("%s_prodVtxYMax%s" ,GetName(),str),"",100,0,10);
454     fhQA[kCutProdVtxZMin] [i] = new TH1F(Form("%s_prodVtxZMin%s" ,GetName(),str),"",100,0,10);
455     fhQA[kCutProdVtxZMax] [i] = new TH1F(Form("%s_prodVtxZMax%s" ,GetName(),str),"",100,0,10);
456     fhQA[kCutDecVtxXMin]  [i] = new TH1F(Form("%s_decVtxXMin%s"  ,GetName(),str),"",100,0,10);
457     fhQA[kCutDecVtxXMax]  [i] = new TH1F(Form("%s_decVtxXMax%s"  ,GetName(),str),"",100,0,10);
458     fhQA[kCutDecVtxYMin]  [i] = new TH1F(Form("%s_decVtxYMin%s"  ,GetName(),str),"",100,0,10);
459     fhQA[kCutDecVtxYMax]  [i] = new TH1F(Form("%s_decVtxYMax%s"  ,GetName(),str),"",100,0,10);
460     fhQA[kCutDecVtxZMin]  [i] = new TH1F(Form("%s_decVtxZMin%s"  ,GetName(),str),"",100,0,10);
461     fhQA[kCutDecVtxZMax]  [i] = new TH1F(Form("%s_decVtxZMax%s"  ,GetName(),str),"",100,0,10);
462     fhQA[kCutDecLgthMin]  [i] = new TH1F(Form("%s_decLengthMin%s",GetName(),str),"",100,0,10);
463     fhQA[kCutDecLgthMax]  [i] = new TH1F(Form("%s_decLengthMax%s",GetName(),str),"",100,0,10);
464     fhQA[kCutDecRxyMin]   [i] = new TH1F(Form("%s_decRxyMin%s"   ,GetName(),str),"",100,0,10);
465     fhQA[kCutDecRxyMax]   [i] = new TH1F(Form("%s_decRxyMax%s"   ,GetName(),str),"",100,0,10);
466     fhQA[kCutDecayChannel][i] = new TH1F(Form("%s_decayChannel%s",GetName(),str),"",2,0,2);
467   }
468   for(Int_t i=0; i<kNCuts; i++) fhQA[i][1]->SetLineColor(color);
469 }
470
471
472 //______________________________
473 Bool_t AliCFParticleGenCuts::IsCharged(AliMCParticle *mcPart) {
474   //
475   //check if particle is charged.
476   //
477   TParticle* part = mcPart->Particle();
478   TParticlePDG* pdgPart = part->GetPDG();
479   if(!pdgPart)return kFALSE;
480   if (pdgPart->Charge() == 0) return kFALSE;
481   return kTRUE;
482 }
483 //______________________________
484 Bool_t AliCFParticleGenCuts::IsPrimary(AliMCParticle *mcPart, AliStack *stack) {
485   //
486   //check if particle is primary (standard definition)
487   //
488   if (!stack->IsPhysicalPrimary(mcPart->Label())) return kFALSE ;
489   return kTRUE;
490 }
491 //______________________________
492 Bool_t AliCFParticleGenCuts::IsPrimaryCharged(AliMCParticle *mcPart, AliStack *stack) {
493   //
494   //check if a charged particle is primary (standard definition)
495   //
496   if (!stack->IsPhysicalPrimary(mcPart->Label()) || !IsCharged(mcPart)) return kFALSE ;
497   return kTRUE;
498 }
499 //______________________________
500 Bool_t AliCFParticleGenCuts::IsA(AliMCParticle *mcPart, Int_t pdg, Bool_t abs) {
501   //
502   //Check on the pdg code of the MC particle. if abs=kTRUE then check on the 
503   //absolute value. By default is set to kFALSE.
504   //
505   TParticle* part = mcPart->Particle();
506   Int_t pdgCode = part->GetPdgCode();
507   if (abs) {
508     pdgCode = TMath::Abs(pdgCode);
509     pdg = TMath::Abs(pdg);
510   }
511   if (pdgCode != pdg ) return kFALSE;
512   return kTRUE;
513 }
514 //______________________________
515 void AliCFParticleGenCuts::SetEvtInfo(TObject* mcEvent) {
516   //
517   // Sets pointer to MC event information (AliMCEvent)
518   //
519
520   if (!mcEvent) {
521     AliError("Pointer to MC Event is null !");
522     return;
523   }
524   
525   TString className(mcEvent->ClassName());
526   if (className.CompareTo("AliMCEvent") != 0) {
527     AliError("argument must point to an AliMCEvent !");
528     return ;
529   }
530   
531   fMCInfo = (AliMCEvent*) mcEvent ;
532 }