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