protection for ProdVtxRange2D added
[u/mrichter/AliRoot.git] / CORRFW / AliCFParticleGenCuts.cxx
index df4d0f8..b66f4e7 100644 (file)
 #include "AliCFParticleGenCuts.h"
 #include "TParticle.h"
 #include "TParticlePDG.h"
-#include "AliMCEventHandler.h"
 #include "AliMCEvent.h"
 #include "TObject.h"
 #include "AliStack.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TBits.h"
+#include "TList.h"
+#include "TArrayF.h"
+#include "TDecayChannel.h"
+#include "AliAODMCParticle.h"
+#include "AliAODEvent.h"
 
 ClassImp(AliCFParticleGenCuts)
 
 //______________________________
 AliCFParticleGenCuts::AliCFParticleGenCuts() : 
   AliCFCutBase(),
+  fIsAODMC(0),
   fMCInfo(0x0),
   fRequireIsCharged(0),
+  fRequireIsNeutral(0),
   fRequireIsPrimary(0),
   fRequireIsSecondary(0),
   fRequirePdgCode(0),
+  fRequireAbsolutePdg(0),
+  fProdVtxRange2D(0),
   fPdgCode(0),
   fProdVtxXMin (-1.e+09),
   fProdVtxYMin (-1.e+09),
@@ -54,24 +65,39 @@ AliCFParticleGenCuts::AliCFParticleGenCuts() :
   fDecayVtxXMax( 1.e+09),
   fDecayVtxYMax( 1.e+09),
   fDecayVtxZMax( 1.e+09),
-  fDecayLengthMin(0),
+  fDecayLengthMin(-1.),
   fDecayLengthMax(1.e+09),
-  fDecayRxyMin(0),
-  fDecayRxyMax(1.e+09)
+  fDecayRxyMin(-1),
+  fDecayRxyMax(1.e+09),
+  fDecayChannel(0x0),
+  fhCutStatistics(0x0),
+  fhCutCorrelation(0x0),
+  fCutValues(new TArrayF(kNCuts)),
+  fBitmap(new TBits(0))
 {
   //
   //ctor
   //
+  for (int i=0; i<kNCuts; i++) 
+    for (int j=0; j<kNStepQA; j++) 
+      fhQA[i][j]=0x0;
+
+  for (int j=0; j<kNStepQA; j++)
+    fhProdVtxXY[j]=0x0;
 }
 
 //______________________________
 AliCFParticleGenCuts::AliCFParticleGenCuts(const Char_t* name, const Char_t* title) : 
   AliCFCutBase(name,title),
+  fIsAODMC(0),
   fMCInfo(0x0),
   fRequireIsCharged(0),
+  fRequireIsNeutral(0),
   fRequireIsPrimary(0),
   fRequireIsSecondary(0),
   fRequirePdgCode(0),
+  fRequireAbsolutePdg(0),
+  fProdVtxRange2D(0),
   fPdgCode(0),
   fProdVtxXMin (-1.e+09),
   fProdVtxYMin (-1.e+09),
@@ -85,24 +111,39 @@ AliCFParticleGenCuts::AliCFParticleGenCuts(const Char_t* name, const Char_t* tit
   fDecayVtxXMax( 1.e+09),
   fDecayVtxYMax( 1.e+09),
   fDecayVtxZMax( 1.e+09),
-  fDecayLengthMin(0),
+  fDecayLengthMin(-1.),
   fDecayLengthMax(1.e+09),
-  fDecayRxyMin(0),
-  fDecayRxyMax(1.e+09)
+  fDecayRxyMin(-1.),
+  fDecayRxyMax(1.e+09),
+  fDecayChannel(0x0),
+  fhCutStatistics(0x0),
+  fhCutCorrelation(0x0),
+  fCutValues(new TArrayF(kNCuts)),
+  fBitmap(new TBits(0))
 {
   //
   //ctor
   //
+  for (int i=0; i<kNCuts; i++) 
+    for (int j=0; j<kNStepQA; j++) 
+      fhQA[i][j]=0x0;
+
+  for (int j=0; j<kNStepQA; j++)
+    fhProdVtxXY[j]=0x0;
 }
 
 //______________________________
 AliCFParticleGenCuts::AliCFParticleGenCuts(const AliCFParticleGenCuts& c) : 
   AliCFCutBase(c),
+  fIsAODMC(c.fIsAODMC),
   fMCInfo(c.fMCInfo),
   fRequireIsCharged(c.fRequireIsCharged),
+  fRequireIsNeutral(c.fRequireIsNeutral),
   fRequireIsPrimary(c.fRequireIsPrimary),
   fRequireIsSecondary(c.fRequireIsSecondary),
   fRequirePdgCode(c.fRequirePdgCode),
+  fRequireAbsolutePdg(c.fRequireAbsolutePdg),
+  fProdVtxRange2D(c.fProdVtxRange2D),
   fPdgCode(c.fPdgCode),
   fProdVtxXMin (c.fProdVtxXMin),
   fProdVtxYMin (c.fProdVtxYMin),
@@ -119,11 +160,22 @@ AliCFParticleGenCuts::AliCFParticleGenCuts(const AliCFParticleGenCuts& c) :
   fDecayLengthMin(c.fDecayLengthMin),
   fDecayLengthMax(c.fDecayLengthMin),
   fDecayRxyMin(c.fDecayLengthMin),
-  fDecayRxyMax(c.fDecayLengthMin)
+  fDecayRxyMax(c.fDecayLengthMin),
+  fDecayChannel(c.fDecayChannel),
+  fhCutStatistics(new TH1F(*c.fhCutStatistics)),
+  fhCutCorrelation(new TH2F(*c.fhCutCorrelation)),
+  fCutValues(new TArrayF(*c.fCutValues)),
+  fBitmap(new TBits(*c.fBitmap))
 {
   //
   //copy ctor
   //
+  for (int i=0; i<kNCuts; i++)
+    for (int j=0; j<kNStepQA; j++)
+      fhQA[i][j]=(TH1F*)c.fhQA[i][j]->Clone();
+
+  for (int j=0; j<kNStepQA; j++)
+    fhProdVtxXY[j]=(TH2F*)c.fhProdVtxXY[j]->Clone();
 }
 
 //______________________________
@@ -134,11 +186,15 @@ AliCFParticleGenCuts& AliCFParticleGenCuts::operator=(const AliCFParticleGenCuts
   //
   if (this != &c) {
     AliCFCutBase::operator=(c) ;
+    fIsAODMC=c.fIsAODMC;
     fMCInfo=c.fMCInfo;
     fRequireIsCharged=c.fRequireIsCharged;
+    fRequireIsNeutral=c.fRequireIsNeutral;
     fRequireIsPrimary=c.fRequireIsPrimary;
     fRequireIsSecondary=c.fRequireIsSecondary;
     fRequirePdgCode=c.fRequirePdgCode;
+    fRequireAbsolutePdg=c.fRequireAbsolutePdg;
+    fProdVtxRange2D=c.fProdVtxRange2D;
     fPdgCode=c.fPdgCode;
     fProdVtxXMin=c.fProdVtxXMin;
     fProdVtxYMin=c.fProdVtxYMin;
@@ -156,6 +212,19 @@ AliCFParticleGenCuts& AliCFParticleGenCuts::operator=(const AliCFParticleGenCuts
     fDecayLengthMax=c.fDecayLengthMax;
     fDecayRxyMin=c.fDecayRxyMin;
     fDecayRxyMax=c.fDecayRxyMax;
+    fDecayChannel=c.fDecayChannel;
+    fCutValues=new TArrayF(*c.fCutValues);
+    fBitmap=new TBits(*c.fBitmap);
+    
+    if (fhCutStatistics)  fhCutStatistics =new TH1F(*c.fhCutStatistics) ;
+    if (fhCutCorrelation) fhCutCorrelation=new TH2F(*c.fhCutCorrelation);
+    
+    for (int i=0; i<kNCuts; i++)
+      for (int j=0; j<kNStepQA; j++)
+       fhQA[i][j]=(TH1F*)c.fhQA[i][j]->Clone();
+
+    for (int j=0; j<kNStepQA; j++)
+       fhProdVtxXY[j]=(TH2F*)c.fhProdVtxXY[j]->Clone();
   }
   return *this ;
 }
@@ -168,120 +237,552 @@ Bool_t AliCFParticleGenCuts::IsSelected(TObject* obj) {
   //
   
   if (!obj) return kFALSE ;
-  TString className(obj->ClassName());
-  if (className.CompareTo("AliMCParticle") != 0) {
-    AliError("argument must point to an AliMCParticle !");
-    return kFALSE ;
+
+  if (!fIsAODMC) SelectionBitMap((AliMCParticle*)   obj);
+  else           SelectionBitMap((AliAODMCParticle*)obj);
+
+  if (fIsQAOn) FillHistograms(obj,0);
+
+  for (UInt_t icut=0; icut<fBitmap->GetNbits();icut++)
+    if (!fBitmap->TestBitNumber(icut)) return kFALSE ; 
+  
+  if (fIsQAOn) FillHistograms(obj,1);
+  return kTRUE;
+}
+            
+//__________________________________________________________________________________
+void AliCFParticleGenCuts::SelectionBitMap(AliMCParticle* mcPart)
+{
+  //
+  // test if the track passes the single cuts
+  // and store the information in a bitmap
+  //
+
+  for (UInt_t i=0; i<kNCuts; i++) {
+    fBitmap->SetBitNumber(i,kFALSE);
+    fCutValues->SetAt((Double32_t)0,i) ;
   }
 
-  AliMCParticle* mcPart = (AliMCParticle*) obj ;
+  // fill the cut array
+  Double32_t partVx=(Double32_t)mcPart->Xv();
+  Double32_t partVy=(Double32_t)mcPart->Yv();
+  Double32_t partVz=(Double32_t)mcPart->Zv();
+
+  // calculate the production vertex ellipse
+  Double32_t prodVtxXYmin = 0.;
+  if (fProdVtxXMin>0 && fProdVtxYMin>0)
+       prodVtxXYmin = partVx*partVx/(fProdVtxXMin*fProdVtxXMin) + partVy*partVy/(fProdVtxYMin*fProdVtxYMin);
+  Double32_t prodVtxXYmax = 0.;
+  if(fProdVtxXMax>0 && fProdVtxYMax>0)
+       prodVtxXYmax = partVx*partVx/(fProdVtxXMax*fProdVtxXMax) + partVy*partVy/(fProdVtxYMax*fProdVtxYMax);
+
+  Double32_t decayVx=0.;
+  Double32_t decayVy=0.;
+  Double32_t decayVz=0.;
+  Double32_t decayL=0.;
+  Double32_t decayRxy=0.;
+
   TParticle* part = mcPart->Particle();
-  AliStack *stack=fMCInfo->MCEvent()->Stack();
+  AliStack*  stack = ((AliMCEvent*)fMCInfo)->Stack();
+  TParticle* daughter=0x0;
+  if ( part->GetNDaughters() > 0 ) {
+    daughter = stack->Particle(part->GetFirstDaughter()) ;
+    decayVx=(Double32_t)daughter->Vx();
+    decayVy=(Double32_t)daughter->Vy();
+    decayVz=(Double32_t)daughter->Vz();
+    decayL = TMath::Sqrt(TMath::Power(partVx-decayVx,2) + 
+                        TMath::Power(partVy-decayVy,2) + 
+                        TMath::Power(partVz-decayVz,2) ) ;
+    decayRxy = TMath::Sqrt(TMath::Power(decayVx,2) + TMath::Power(decayVy,2) ) ;
+  }
 
-  // is this particle charged?
-  if ( fRequireIsCharged ) {
-    if(!IsCharged(mcPart))return kFALSE;
+  fCutValues->SetAt(partVx  ,kCutProdVtxXMin) ;
+  fCutValues->SetAt(partVy  ,kCutProdVtxYMin) ;
+  fCutValues->SetAt(partVz  ,kCutProdVtxZMin) ;
+  fCutValues->SetAt(partVx  ,kCutProdVtxXMax) ;
+  fCutValues->SetAt(partVy  ,kCutProdVtxYMax) ;
+  fCutValues->SetAt(partVz  ,kCutProdVtxZMax) ;
+  fCutValues->SetAt(decayVx ,kCutDecVtxXMin)  ;
+  fCutValues->SetAt(decayVy ,kCutDecVtxYMin)  ;
+  fCutValues->SetAt(decayVz ,kCutDecVtxZMin)  ;
+  fCutValues->SetAt(decayVx ,kCutDecVtxXMax)  ;
+  fCutValues->SetAt(decayVy ,kCutDecVtxYMax)  ;
+  fCutValues->SetAt(decayVz ,kCutDecVtxZMax)  ;
+  fCutValues->SetAt(decayL  ,kCutDecLgthMin)  ;
+  fCutValues->SetAt(decayL  ,kCutDecLgthMax)  ;
+  fCutValues->SetAt(decayRxy,kCutDecRxyMin)   ;
+  fCutValues->SetAt(decayRxy,kCutDecRxyMax)   ;
+  
+  // cut on charge
+  if ( fRequireIsCharged || fRequireIsNeutral ) {
+    if (fRequireIsCharged &&  IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
+    if (fRequireIsNeutral && !IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
   } 
+  else fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
   
-  // primary cuts
-  if ( fRequireIsPrimary ) {
-     if(!IsPrimary(mcPart,stack))return kFALSE;
+  // cut on primary/secondary
+  if ( fRequireIsPrimary || fRequireIsSecondary) {
+    if (fRequireIsPrimary   &&  IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
+    if (fRequireIsSecondary && !IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
   } 
+  else fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
+  
+  // cut on PDG code
+  if ( fRequirePdgCode ) {
+    if (IsA(mcPart,fPdgCode,fRequireAbsolutePdg)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
+  }
+  else fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
+  
+  // cut on decay channel
+  if ( fDecayChannel ) {
+    Bool_t goodDecay = kTRUE ;
+    Short_t nDaughters = mcPart->Particle()->GetNDaughters() ;
+    if (nDaughters != fDecayChannel->NDaughters()) goodDecay = kFALSE ;
+    //now number of daughters is OK
+    if (goodDecay) {
+      // now check if decay channel is respected
+      // first try
+      for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) {
+       TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(iDaughter)) ;
+       if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
+      }
+      if (!goodDecay) {
+       //second try inverting the order of the daughters
+       goodDecay = kTRUE ;
+       for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) {
+         TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(nDaughters-(iDaughter+1))) ;
+         if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
+       }
+      }
+      if (!goodDecay && fRequireAbsolutePdg) {
+       //now tries inverting the sign of the daughters in case the anti-particle is also looked at
+       // third try
+       goodDecay = kTRUE ;
+       for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) {
+         TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(iDaughter)) ;
+         if (daug->GetPdgCode() != -fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
+       }
+       if (!goodDecay) {
+         //fourth try inverting the order of the daughters
+         goodDecay = kTRUE ;
+         for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) {
+           TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(nDaughters-(iDaughter+1))) ;
+           if (daug->GetPdgCode() != -fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
+         }
+       }
+      } //end check anti-particle
+    } //end # daughters OK
+    fCutValues->SetAt((Double32_t)goodDecay,kCutDecayChannel) ;
+  } //end require decay channel
+  else fCutValues->SetAt((Double32_t)kTRUE,kCutDecayChannel);
+  
+  
+  // now array of cut is build, fill the bitmap consequently
+  Int_t iCutBit = -1;
+  if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
+
+  ++iCutBit;
+  if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxXMin)
+    || ( fProdVtxRange2D && (fProdVtxXMin>0 && fProdVtxYMin>0) && prodVtxXYmin >= 1)
+    || ( fProdVtxRange2D && (fProdVtxXMin<=0 || fProdVtxYMin<=0) ) )
+   fBitmap->SetBitNumber(iCutBit,kTRUE);
+
+  ++iCutBit;
+  if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxXMax)
+    || ( fProdVtxRange2D && (fProdVtxXMax>0 && fProdVtxYMax>0) && prodVtxXYmax <= 1)
+    || ( fProdVtxRange2D && (fProdVtxXMax<=0 || fProdVtxYMax<=0) ) )
+  fBitmap->SetBitNumber(iCutBit,kTRUE);
+
+  ++iCutBit;
+  if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxYMin)
+    || ( fProdVtxRange2D &&  (fProdVtxXMin>0 && fProdVtxYMin>0) && prodVtxXYmin >= 1)
+    || ( fProdVtxRange2D &&  (fProdVtxXMin<=0 || fProdVtxYMin<=0) ) )
+  fBitmap->SetBitNumber(iCutBit,kTRUE);
 
-  //secondary cut
-  if ( fRequireIsSecondary && part->IsPrimary() ) return kFALSE ;
+  ++iCutBit;
+  if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxYMax)
+    || ( fProdVtxRange2D && (fProdVtxXMax>0 && fProdVtxYMax>0) && prodVtxXYmax <= 1)
+    || ( fProdVtxRange2D && (fProdVtxXMax<=0 || fProdVtxYMax<=0) ) )
+  fBitmap->SetBitNumber(iCutBit,kTRUE);
+
+  if ( fCutValues->At(++iCutBit) > fProdVtxZMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) < fProdVtxZMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) > fDecayVtxXMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) < fDecayVtxXMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) > fDecayVtxYMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) < fDecayVtxYMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) > fDecayVtxZMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) < fDecayVtxZMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) > fDecayLengthMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) < fDecayLengthMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) > fDecayRxyMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) < fDecayRxyMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) != 0 )             fBitmap->SetBitNumber(iCutBit,kTRUE);
+}
+
+//__________________________________________________________________________________
+void AliCFParticleGenCuts::SelectionBitMap(AliAODMCParticle* mcPart)
+{
+  //
+  // test if the track passes the single cuts
+  // and store the information in a bitmap
+  //
   
-  //PDG code cut
-  if ( fRequirePdgCode){
-    if(!IsA(mcPart,fPdgCode))  return kFALSE ;
+  for (UInt_t i=0; i<kNCuts; i++) {
+    fBitmap->SetBitNumber(i,kFALSE);
+    fCutValues->SetAt((Double32_t)0,i) ;
   }
-  // production vertex cuts
-  Double32_t partVx=(Double32_t)part->Vx();
-  Double32_t partVy=(Double32_t)part->Vy();
-  Double32_t partVz=(Double32_t)part->Vz();
-  if ( partVx < fProdVtxXMin || partVx > fProdVtxXMax ) return kFALSE ;
-  if ( partVy < fProdVtxYMin || partVy > fProdVtxYMax ) return kFALSE ;
-  if ( partVz < fProdVtxZMin || partVz > fProdVtxZMax ) return kFALSE ;
-
-  //decay vertex cuts
-  if ( part->GetNDaughters() > 0 ) {
-    TParticle* daughter = fMCInfo->MCEvent()->Stack()->Particle(part->GetFirstDaughter()) ;
-    Double32_t decayVx=(Double32_t)daughter->Vx();
-    Double32_t decayVy=(Double32_t)daughter->Vy();
-    Double32_t decayVz=(Double32_t)daughter->Vz();
-    if ( decayVx < fDecayVtxXMin || decayVx > fDecayVtxXMax ) return kFALSE ;
-    if ( decayVy < fDecayVtxYMin || decayVy > fDecayVtxYMax ) return kFALSE ;
-    if ( decayVz < fDecayVtxZMin || decayVz > fDecayVtxZMax ) return kFALSE ;
-
-    //decay length cut
-    Double32_t decayL = TMath::Sqrt(TMath::Power(partVx-decayVx,2) + 
-                                   TMath::Power(partVy-decayVy,2) + 
-                                   TMath::Power(partVz-decayVz,2) ) ;
-    if (decayL < fDecayLengthMin || decayL > fDecayLengthMax) return kFALSE ;
-
-    Double32_t decayRxy = TMath::Sqrt(TMath::Power(decayVx,2) + 
-                                     TMath::Power(decayVy,2) ) ;
-    if (decayRxy < fDecayRxyMin || decayRxy > fDecayRxyMax) return kFALSE ;
+
+  // fill the cut array
+  Double32_t partVx=(Double32_t)mcPart->Xv();
+  Double32_t partVy=(Double32_t)mcPart->Yv();
+  Double32_t partVz=(Double32_t)mcPart->Zv();
+
+  // calculate the production vertex ellipse
+  Double32_t prodVtxXYmin = 0.;
+  if (fProdVtxXMin!=0 && fProdVtxYMin!=0)
+       prodVtxXYmin = partVx*partVx/(fProdVtxXMin*fProdVtxXMin) + partVy*partVy/(fProdVtxYMin*fProdVtxYMin);
+  Double32_t prodVtxXYmax = 0.;
+  if(fProdVtxXMax!=0 && fProdVtxYMax!=0)
+       prodVtxXYmax = partVx*partVx/(fProdVtxXMax*fProdVtxXMax) + partVy*partVy/(fProdVtxYMax*fProdVtxYMax);
+
+  Double32_t decayVx=0.;
+  Double32_t decayVy=0.;
+  Double32_t decayVz=0.;
+  Double32_t decayL=0.;
+  Double32_t decayRxy=0.;
+
+  AliAODEvent* aod = dynamic_cast<AliAODEvent*>(fMCInfo);
+  TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
+  AliAODMCParticle* daughter = 0x0 ;
+
+  if (mcPart->GetDaughter(0)>0) {
+    daughter = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(0)));
+
+    decayVx=(Double32_t)daughter->Xv();
+    decayVy=(Double32_t)daughter->Yv();
+    decayVz=(Double32_t)daughter->Zv();
+    decayL = TMath::Sqrt(TMath::Power(partVx-decayVx,2) + 
+                        TMath::Power(partVy-decayVy,2) + 
+                        TMath::Power(partVz-decayVz,2) ) ;
+    decayRxy = TMath::Sqrt(TMath::Power(decayVx,2) + TMath::Power(decayVy,2) ) ;
+    
+  }
+  
+  fCutValues->SetAt(partVx  ,kCutProdVtxXMin) ;
+  fCutValues->SetAt(partVy  ,kCutProdVtxYMin) ;
+  fCutValues->SetAt(partVz  ,kCutProdVtxZMin) ;
+  fCutValues->SetAt(partVx  ,kCutProdVtxXMax) ;
+  fCutValues->SetAt(partVy  ,kCutProdVtxYMax) ;
+  fCutValues->SetAt(partVz  ,kCutProdVtxZMax) ;
+  fCutValues->SetAt(decayVx ,kCutDecVtxXMin)  ;
+  fCutValues->SetAt(decayVy ,kCutDecVtxYMin)  ;
+  fCutValues->SetAt(decayVz ,kCutDecVtxZMin)  ;
+  fCutValues->SetAt(decayVx ,kCutDecVtxXMax)  ;
+  fCutValues->SetAt(decayVy ,kCutDecVtxYMax)  ;
+  fCutValues->SetAt(decayVz ,kCutDecVtxZMax)  ;
+  fCutValues->SetAt(decayL  ,kCutDecLgthMin)  ;
+  fCutValues->SetAt(decayL  ,kCutDecLgthMax)  ;
+  fCutValues->SetAt(decayRxy,kCutDecRxyMin)   ;
+  fCutValues->SetAt(decayRxy,kCutDecRxyMax)   ;
+  
+  // cut on charge
+  if ( fRequireIsCharged || fRequireIsNeutral ) {
+    if (fRequireIsCharged &&  IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
+    if (fRequireIsNeutral && !IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
+  } 
+  else fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
+  
+  // cut on primary/secondary
+  if ( fRequireIsPrimary || fRequireIsSecondary) {
+    if (fRequireIsPrimary   &&  IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
+    if (fRequireIsSecondary && !IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
+  } 
+  else fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
+  
+  // cut on PDG code
+  if ( fRequirePdgCode ) {
+    if (IsA(mcPart,fPdgCode,fRequireAbsolutePdg)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
+  }
+  else fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
+  
+  // cut on decay channel
+  if ( fDecayChannel ) {
+    Bool_t goodDecay = kTRUE ;
+    Short_t nDaughters = 0 ;
+    if (mcPart->GetDaughter(0) >=0) nDaughters = 1 + mcPart->GetDaughter(1) - mcPart->GetDaughter(0);
+    
+    if (nDaughters != fDecayChannel->NDaughters()) goodDecay = kFALSE ;
+    if (goodDecay) {
+      // now check if decay channel is respected
+      // first try
+      for (Int_t iDaughter = 0 ; iDaughter<nDaughters; iDaughter++) {
+       AliAODMCParticle* daug = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(0)+iDaughter));
+       if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
+      }
+      if (!goodDecay) {
+       //second try inverting the order of the daughters
+       goodDecay = kTRUE ;
+       for (Int_t iDaughter = 0 ; iDaughter<nDaughters; iDaughter++) {
+         AliAODMCParticle* daug = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(1)-iDaughter));
+         if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
+       }
+      }
+    }
+    fCutValues->SetAt((Double32_t)goodDecay,kCutDecayChannel) ;
   }
+  else fCutValues->SetAt((Double32_t)kTRUE,kCutDecayChannel);
+  
+  
+  // now array of cut is build, fill the bitmap consequently
+  Int_t iCutBit = -1;
+  if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
+
+  ++iCutBit;
+  if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxXMin)
+    || ( fProdVtxRange2D && prodVtxXYmin >= 1))     fBitmap->SetBitNumber(iCutBit,kTRUE);
 
+  ++iCutBit;
+  if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxXMax)
+    || ( fProdVtxRange2D && prodVtxXYmax <= 1))     fBitmap->SetBitNumber(iCutBit,kTRUE);
 
-  return kTRUE ;
+  ++iCutBit;
+  if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxYMin)
+    || ( fProdVtxRange2D && prodVtxXYmin >= 1))     fBitmap->SetBitNumber(iCutBit,kTRUE);
+
+  ++iCutBit;
+  if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxYMax)
+    || ( fProdVtxRange2D && prodVtxXYmax <= 1))     fBitmap->SetBitNumber(iCutBit,kTRUE);
+
+  if ( fCutValues->At(++iCutBit) > fProdVtxZMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) < fProdVtxZMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) > fDecayVtxXMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) < fDecayVtxXMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) > fDecayVtxYMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) < fDecayVtxYMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) > fDecayVtxZMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) < fDecayVtxZMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) > fDecayLengthMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) < fDecayLengthMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) > fDecayRxyMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) < fDecayRxyMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
+  if ( fCutValues->At(++iCutBit) != 0 )             fBitmap->SetBitNumber(iCutBit,kTRUE);
 }
+
+
+//__________________________________________________________________________________
+void AliCFParticleGenCuts::FillHistograms(TObject* /*obj*/, Bool_t afterCuts)
+{
+  //
+  // fill the QA histograms
+  //
+
+  for (int iCutNumber = 0; iCutNumber < kNCuts; iCutNumber++) 
+    fhQA[iCutNumber][afterCuts]->Fill(fCutValues->At(iCutNumber));
+
+  fhProdVtxXY[afterCuts]->Fill(fCutValues->At(4),fCutValues->At(5));
+
+  // fill cut statistics and cut correlation histograms with information from the bitmap
+  if (afterCuts) return;
+
+  // Number of single cuts in this class
+  UInt_t ncuts = fBitmap->GetNbits();
+  for(UInt_t bit=0; bit<ncuts;bit++) {
+    if (!fBitmap->TestBitNumber(bit)) {
+      fhCutStatistics->Fill(bit+1);
+      for (UInt_t bit2=bit; bit2<ncuts;bit2++) {
+       if (!fBitmap->TestBitNumber(bit2)) 
+         fhCutCorrelation->Fill(bit+1,bit2+1);
+      }
+    }
+  }
+}
+
+//__________________________________________________________________________________
+void AliCFParticleGenCuts::AddQAHistograms(TList *qaList) {
+  //
+  // saves the histograms in a TList
+  //
+
+  DefineHistograms();
+
+  qaList->Add(fhCutStatistics);
+  qaList->Add(fhCutCorrelation);
+
+  for (Int_t j=0; j<kNStepQA; j++) {
+    qaList->Add(fhProdVtxXY[j]);
+    for(Int_t i=0; i<kNCuts; i++)
+      qaList->Add(fhQA[i][j]);
+  }
+}
+
+//__________________________________________________________________________________
+void AliCFParticleGenCuts::DefineHistograms() {
+  //
+  // histograms for cut variables, cut statistics and cut correlations
+  //
+  Int_t color = 2;
+
+  // book cut statistics and cut correlation histograms
+  fhCutStatistics = new TH1F(Form("%s_cut_statistics",GetName()),"",kNCuts,0.5,kNCuts+0.5);
+  fhCutStatistics->SetLineWidth(2);
+  int k = 1;
+  fhCutStatistics->GetXaxis()->SetBinLabel(k,"charge")     ; k++;
+  fhCutStatistics->GetXaxis()->SetBinLabel(k,"prim/sec")   ; k++;
+  fhCutStatistics->GetXaxis()->SetBinLabel(k,"PDG")        ; k++;
+  fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxXMin")    ; k++;
+  fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxXMax")    ; k++;
+  fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxYMin")    ; k++;
+  fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxYMax")    ; k++;
+  fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxZMin")    ; k++;
+  fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMax")    ; k++;
+  fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecXMin")    ; k++;
+  fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecXMax")    ; k++;
+  fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecYMin")    ; k++;
+  fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecYMax")    ; k++;
+  fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMin")    ; k++;
+  fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMax")    ; k++;
+  fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecLgthMin") ; k++;
+  fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecLgthMax") ; k++;
+  fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMin")  ; k++;
+  fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMax")  ; k++;
+  fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecChannel") ; k++;
+
+
+  fhCutCorrelation = new TH2F(Form("%s_cut_correlation",GetName()),"",kNCuts,0.5,kNCuts+0.5,kNCuts,0.5,kNCuts+0.5);
+  fhCutCorrelation->SetLineWidth(2);
+  for (k=1; k<=kNCuts; k++) {
+    fhCutCorrelation->GetXaxis()->SetBinLabel(k,fhCutStatistics->GetXaxis()->GetBinLabel(k));
+    fhCutCorrelation->GetYaxis()->SetBinLabel(k,fhCutStatistics->GetXaxis()->GetBinLabel(k));
+  }
+
+  Char_t str[256];
+  for (int i=0; i<kNStepQA; i++) {
+    if (i==0) sprintf(str," ");
+    else sprintf(str,"_cut");
+    fhQA[kCutCharge]      [i] = new TH1F(Form("%s_charge%s"      ,GetName(),str),"",2,0,2);
+    fhQA[kCutPrimSec]     [i] = new TH1F(Form("%s_primSec%s"     ,GetName(),str),"",2,0,2);
+    fhQA[kCutPDGCode]     [i] = new TH1F(Form("%s_pdgCode%s"     ,GetName(),str),"",2,0,2);
+    fhQA[kCutProdVtxXMin] [i] = new TH1F(Form("%s_prodVtxXMin%s" ,GetName(),str),"",100,-10,10);
+    fhQA[kCutProdVtxXMax] [i] = new TH1F(Form("%s_prodVtxXMax%s" ,GetName(),str),"",100,-10,10);
+    fhQA[kCutProdVtxYMin] [i] = new TH1F(Form("%s_prodVtxYMin%s" ,GetName(),str),"",100,-10,10);
+    fhQA[kCutProdVtxYMax] [i] = new TH1F(Form("%s_prodVtxYMax%s" ,GetName(),str),"",100,-10,10);
+    fhQA[kCutProdVtxZMin] [i] = new TH1F(Form("%s_prodVtxZMin%s" ,GetName(),str),"",100,-10,10);
+    fhQA[kCutProdVtxZMax] [i] = new TH1F(Form("%s_prodVtxZMax%s" ,GetName(),str),"",100,-10,10);
+    fhQA[kCutDecVtxXMin]  [i] = new TH1F(Form("%s_decVtxXMin%s"  ,GetName(),str),"",100,0,10);
+    fhQA[kCutDecVtxXMax]  [i] = new TH1F(Form("%s_decVtxXMax%s"  ,GetName(),str),"",100,0,10);
+    fhQA[kCutDecVtxYMin]  [i] = new TH1F(Form("%s_decVtxYMin%s"  ,GetName(),str),"",100,0,10);
+    fhQA[kCutDecVtxYMax]  [i] = new TH1F(Form("%s_decVtxYMax%s"  ,GetName(),str),"",100,0,10);
+    fhQA[kCutDecVtxZMin]  [i] = new TH1F(Form("%s_decVtxZMin%s"  ,GetName(),str),"",100,0,10);
+    fhQA[kCutDecVtxZMax]  [i] = new TH1F(Form("%s_decVtxZMax%s"  ,GetName(),str),"",100,0,10);
+    fhQA[kCutDecLgthMin]  [i] = new TH1F(Form("%s_decLengthMin%s",GetName(),str),"",100,0,10);
+    fhQA[kCutDecLgthMax]  [i] = new TH1F(Form("%s_decLengthMax%s",GetName(),str),"",100,0,10);
+    fhQA[kCutDecRxyMin]   [i] = new TH1F(Form("%s_decRxyMin%s"   ,GetName(),str),"",100,0,10);
+    fhQA[kCutDecRxyMax]   [i] = new TH1F(Form("%s_decRxyMax%s"   ,GetName(),str),"",100,0,10);
+    fhQA[kCutDecayChannel][i] = new TH1F(Form("%s_decayChannel%s",GetName(),str),"",2,0,2);
+    fhProdVtxXY                  [i] = new TH2F(Form("%s_prodVtxXY%s"   ,GetName(),str),"",100,0,10,100,0,10);
+    fhProdVtxXY                  [i] ->GetXaxis()->SetTitle("x_{production vertex}");
+    fhProdVtxXY                  [i] ->GetYaxis()->SetTitle("y_{production vertex}");
+    fhQA[kCutProdVtxXMax] [i] ->GetXaxis()->SetTitle("x_{production vertex}");
+    fhQA[kCutProdVtxYMax] [i] ->GetXaxis()->SetTitle("y_{production vertex}");
+  }
+  for(Int_t i=0; i<kNCuts; i++) fhQA[i][1]->SetLineColor(color);
+}
+
+
 //______________________________
-Bool_t AliCFParticleGenCuts::IsCharged(AliMCParticle *mcPart) {
+Bool_t AliCFParticleGenCuts::IsCharged(AliVParticle *mcPart) {
   //
   //check if particle is charged.
   //
-  TParticle* part = mcPart->Particle();
-  TParticlePDG* pdgPart = part->GetPDG();
-  if(!pdgPart)return kFALSE;
-  if (pdgPart->Charge() == 0) return kFALSE;
+  if (mcPart->Charge()==0) return kFALSE;
+  return kTRUE;
+}
+//______________________________
+Bool_t AliCFParticleGenCuts::IsPrimary(AliMCParticle *mcPart) {
+  //
+  //check if particle is primary (standard definition)
+  //
+  
+  AliStack* stack = ((AliMCEvent*)fMCInfo)->Stack();
+
+  if (!stack->IsPhysicalPrimary(mcPart->GetLabel())) return kFALSE;
   return kTRUE;
 }
 //______________________________
-Bool_t AliCFParticleGenCuts::IsPrimary(AliMCParticle *mcPart, AliStack *stack) {
+Bool_t AliCFParticleGenCuts::IsPrimary(AliAODMCParticle *mcPart) {
   //
   //check if particle is primary (standard definition)
   //
-  if (!stack->IsPhysicalPrimary(mcPart->Label())) return kFALSE ;
+  
+  if (!mcPart->IsPhysicalPrimary()) return kFALSE;
   return kTRUE;
 }
 //______________________________
-Bool_t AliCFParticleGenCuts::IsPrimaryCharged(AliMCParticle *mcPart, AliStack *stack) {
+Bool_t AliCFParticleGenCuts::IsPrimaryCharged(AliVParticle *mcPart) {
   //
   //check if a charged particle is primary (standard definition)
   //
-  if (!stack->IsPhysicalPrimary(mcPart->Label()) || !IsCharged(mcPart)) return kFALSE ;
+
+  if (!fIsAODMC) {
+    if (!IsPrimary((AliMCParticle*)mcPart) || !IsCharged(mcPart)) return kFALSE ;
+  }
+  else {
+    if (!IsPrimary((AliAODMCParticle*)mcPart) || !IsCharged(mcPart)) return kFALSE ;
+  }
   return kTRUE;
 }
 //______________________________
 Bool_t AliCFParticleGenCuts::IsA(AliMCParticle *mcPart, Int_t pdg, Bool_t abs) {
   //
   //Check on the pdg code of the MC particle. if abs=kTRUE then check on the 
-  //absolute value. By default is set to kFALSE.
+  //absolute value. 
   //
   TParticle* part = mcPart->Particle();
   Int_t pdgCode = part->GetPdgCode();
-  if(abs)pdgCode=TMath::Abs(pdgCode);
-  if(pdgCode != pdg )return kFALSE;
+
+  if (abs) {
+    pdgCode = TMath::Abs(pdgCode);
+    pdg = TMath::Abs(pdg);
+  }
+  if (pdgCode != pdg ) return kFALSE;
   return kTRUE;
 }
 //______________________________
-void AliCFParticleGenCuts::SetEvtInfo(TObject* mcInfo) {
+Bool_t AliCFParticleGenCuts::IsA(AliAODMCParticle *mcPart, Int_t pdg, Bool_t abs) {
   //
-  // Sets pointer to MC event information (AliMCEventHandler)
+  //Check on the pdg code of the MC particle. if abs=kTRUE then check on the 
+  //absolute value. 
+  //
+  Int_t pdgCode = mcPart->GetPdgCode();
+  
+  if (abs) {
+    pdgCode = TMath::Abs(pdgCode);
+    pdg = TMath::Abs(pdg);
+  }
+  if (pdgCode != pdg ) return kFALSE;
+  return kTRUE;
+}
+//______________________________
+void AliCFParticleGenCuts::SetMCEventInfo(const TObject* mcEvent) {
+  //
+  // Sets pointer to MC event information (AliMCEvent)
   //
 
-  if (!mcInfo) {
-    AliError("Pointer to MC Event Handler is null !");
+  if (!mcEvent) {
+    AliError("Pointer to MC Event is null !");
     return;
   }
   
-  TString className(mcInfo->ClassName());
-  if (className.CompareTo("AliMCEventHandler") != 0) {
-    AliError("argument must point to an AliMCEventHandler !");
+  TString className(mcEvent->ClassName());
+  if (className.CompareTo("AliMCEvent") != 0 && className.CompareTo("AliAODEvent") != 0) {
+    AliError("argument must point to an AliMCEvent or an AliAODEvent !");
     return ;
   }
-  
-  fMCInfo = (AliMCEventHandler*) mcInfo ;
+
+  fMCInfo = (AliVEvent*)mcEvent ;
 }