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 }