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 }