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   fPdgCode(0),
  54   fProdVtxXMin (-1.e+09),
  55   fProdVtxYMin (-1.e+09),
  56   fProdVtxZMin (-1.e+09),
  57   fProdVtxXMax ( 1.e+09),
  58   fProdVtxYMax ( 1.e+09),
  59   fProdVtxZMax ( 1.e+09),
  60   fDecayVtxXMin(-1.e+09),
  61   fDecayVtxYMin(-1.e+09),
  62   fDecayVtxZMin(-1.e+09),
  63   fDecayVtxXMax( 1.e+09),
  64   fDecayVtxYMax( 1.e+09),
  65   fDecayVtxZMax( 1.e+09),
  66   fDecayLengthMin(-1.),
  67   fDecayLengthMax(1.e+09),
  68   fDecayRxyMin(-1),
  69   fDecayRxyMax(1.e+09),
  70   fDecayChannel(0x0),
  71   fhCutStatistics(0x0),
  72   fhCutCorrelation(0x0),
  73   fCutValues(new TArrayF(kNCuts)),
  74   fBitmap(new TBits(0))
  75 {
  76   //
  77   //ctor
  78   //
  79   for (int i=0; i<kNCuts; i++)
  80     for (int j=0; j<kNStepQA; j++)
  81       fhQA[i][j]=0x0;
  82 }
  83
  84 //______________________________
  85 AliCFParticleGenCuts::AliCFParticleGenCuts(const Char_t* name, const Char_t* title) :
  86   AliCFCutBase(name,title),
  87   fIsAODMC(0),
  88   fMCInfo(0x0),
  89   fRequireIsCharged(0),
  90   fRequireIsNeutral(0),
  91   fRequireIsPrimary(0),
  92   fRequireIsSecondary(0),
  93   fRequirePdgCode(0),
  94   fPdgCode(0),
  95   fProdVtxXMin (-1.e+09),
  96   fProdVtxYMin (-1.e+09),
  97   fProdVtxZMin (-1.e+09),
  98   fProdVtxXMax ( 1.e+09),
  99   fProdVtxYMax ( 1.e+09),
 100   fProdVtxZMax ( 1.e+09),
 101   fDecayVtxXMin(-1.e+09),
 102   fDecayVtxYMin(-1.e+09),
 103   fDecayVtxZMin(-1.e+09),
 104   fDecayVtxXMax( 1.e+09),
 105   fDecayVtxYMax( 1.e+09),
 106   fDecayVtxZMax( 1.e+09),
 107   fDecayLengthMin(-1.),
 108   fDecayLengthMax(1.e+09),
 109   fDecayRxyMin(-1.),
 110   fDecayRxyMax(1.e+09),
 111   fDecayChannel(0x0),
 112   fhCutStatistics(0x0),
 113   fhCutCorrelation(0x0),
 114   fCutValues(new TArrayF(kNCuts)),
 115   fBitmap(new TBits(0))
 116 {
 117   //
 118   //ctor
 119   //
 120   for (int i=0; i<kNCuts; i++)
 121     for (int j=0; j<kNStepQA; j++)
 122       fhQA[i][j]=0x0;
 123 }
 124
 125 //______________________________
 126 AliCFParticleGenCuts::AliCFParticleGenCuts(const AliCFParticleGenCuts& c) :
 127   AliCFCutBase(c),
 128   fIsAODMC(c.fIsAODMC),
 129   fMCInfo(c.fMCInfo),
 130   fRequireIsCharged(c.fRequireIsCharged),
 131   fRequireIsNeutral(c.fRequireIsNeutral),
 132   fRequireIsPrimary(c.fRequireIsPrimary),
 133   fRequireIsSecondary(c.fRequireIsSecondary),
 134   fRequirePdgCode(c.fRequirePdgCode),
 135   fPdgCode(c.fPdgCode),
 136   fProdVtxXMin (c.fProdVtxXMin),
 137   fProdVtxYMin (c.fProdVtxYMin),
 138   fProdVtxZMin (c.fProdVtxZMin),
 139   fProdVtxXMax (c.fProdVtxXMax),
 140   fProdVtxYMax (c.fProdVtxYMax),
 141   fProdVtxZMax (c.fProdVtxZMax),
 142   fDecayVtxXMin(c.fDecayVtxXMin),
 143   fDecayVtxYMin(c.fDecayVtxYMin),
 144   fDecayVtxZMin(c.fDecayVtxZMin),
 145   fDecayVtxXMax(c.fDecayVtxXMax),
 146   fDecayVtxYMax(c.fDecayVtxYMax),
 147   fDecayVtxZMax(c.fDecayVtxZMax),
 148   fDecayLengthMin(c.fDecayLengthMin),
 149   fDecayLengthMax(c.fDecayLengthMin),
 150   fDecayRxyMin(c.fDecayLengthMin),
 151   fDecayRxyMax(c.fDecayLengthMin),
 152   fDecayChannel(c.fDecayChannel),
 153   fhCutStatistics(new TH1F(*c.fhCutStatistics)),
 154   fhCutCorrelation(new TH2F(*c.fhCutCorrelation)),
 155   fCutValues(new TArrayF(*c.fCutValues)),
 156   fBitmap(new TBits(*c.fBitmap))
 157 {
 158   //
 159   //copy ctor
 160   //
 161   for (int i=0; i<kNCuts; i++)
 162     for (int j=0; j<kNStepQA; j++)
 163       fhQA[i][j]=(TH1F*)c.fhQA[i][j]->Clone();
 164 }
 165
 166 //______________________________
 167 AliCFParticleGenCuts& AliCFParticleGenCuts::operator=(const AliCFParticleGenCuts& c)
 168 {
 169   //
 170   // Assignment operator
 171   //
 172   if (this != &c) {
 173     AliCFCutBase::operator=(c) ;
 174     fIsAODMC=c.fIsAODMC;
 175     fMCInfo=c.fMCInfo;
 176     fRequireIsCharged=c.fRequireIsCharged;
 177     fRequireIsNeutral=c.fRequireIsNeutral;
 178     fRequireIsPrimary=c.fRequireIsPrimary;
 179     fRequireIsSecondary=c.fRequireIsSecondary;
 180     fRequirePdgCode=c.fRequirePdgCode;
 181     fPdgCode=c.fPdgCode;
 182     fProdVtxXMin=c.fProdVtxXMin;
 183     fProdVtxYMin=c.fProdVtxYMin;
 184     fProdVtxZMin=c.fProdVtxZMin;
 185     fProdVtxXMax=c.fProdVtxXMax;
 186     fProdVtxYMax=c.fProdVtxYMax;
 187     fProdVtxZMax=c.fProdVtxZMax;
 188     fDecayVtxXMin=c.fDecayVtxXMin;
 189     fDecayVtxYMin=c.fDecayVtxYMin;
 190     fDecayVtxZMin=c.fDecayVtxZMin;
 191     fDecayVtxXMax=c.fDecayVtxXMax;
 192     fDecayVtxYMax=c.fDecayVtxYMax;
 193     fDecayVtxZMax=c.fDecayVtxZMax;
 194     fDecayLengthMin=c.fDecayVtxZMax;
 195     fDecayLengthMax=c.fDecayLengthMax;
 196     fDecayRxyMin=c.fDecayRxyMin;
 197     fDecayRxyMax=c.fDecayRxyMax;
 198     fDecayChannel=c.fDecayChannel;
 199     fCutValues=new TArrayF(*c.fCutValues);
 200     fBitmap=new TBits(*c.fBitmap);
 201
 202     if (fhCutStatistics)  fhCutStatistics =new TH1F(*c.fhCutStatistics) ;
 203     if (fhCutCorrelation) fhCutCorrelation=new TH2F(*c.fhCutCorrelation);
 204
 205     for (int i=0; i<kNCuts; i++)
 206       for (int j=0; j<kNStepQA; j++)
 207         fhQA[i][j]=(TH1F*)c.fhQA[i][j]->Clone();
 208   }
 209   return *this ;
 210 }
 211
 212 //______________________________
 213 Bool_t AliCFParticleGenCuts::IsSelected(TObject* obj) {
 214   //
 215   // check if selections on 'obj' are passed
 216   // 'obj' must be an AliMCParticle
 217   //
 218
 219   if (!obj) return kFALSE ;
 220
 221   if (!fIsAODMC) SelectionBitMap((AliMCParticle*)   obj);
 222   else           SelectionBitMap((AliAODMCParticle*)obj);
 223
 224   if (fIsQAOn) FillHistograms(obj,0);
 225
 226   for (UInt_t icut=0; icut<fBitmap->GetNbits();icut++)
 227     if (!fBitmap->TestBitNumber(icut)) return kFALSE ;
 228
 229   if (fIsQAOn) FillHistograms(obj,1);
 230   return kTRUE;
 231 }
 232
 233 //__________________________________________________________________________________
 234 void AliCFParticleGenCuts::SelectionBitMap(AliMCParticle* mcPart)
 235 {
 236   //
 237   // test if the track passes the single cuts
 238   // and store the information in a bitmap
 239   //
 240
 241   for (UInt_t i=0; i<kNCuts; i++) {
 242     fBitmap->SetBitNumber(i,kFALSE);
 243     fCutValues->SetAt((Double32_t)0,i) ;
 244   }
 245
 246   // fill the cut array
 247   Double32_t partVx=(Double32_t)mcPart->Xv();
 248   Double32_t partVy=(Double32_t)mcPart->Yv();
 249   Double32_t partVz=(Double32_t)mcPart->Zv();
 250
 251   Double32_t decayVx=0.;
 252   Double32_t decayVy=0.;
 253   Double32_t decayVz=0.;
 254   Double32_t decayL=0.;
 255   Double32_t decayRxy=0.;
 256
 257   TParticle* part = mcPart->Particle();
 258   AliStack*  stack = ((AliMCEvent*)fMCInfo)->Stack();
 259   TParticle* daughter=0x0;
 260   if ( part->GetNDaughters() > 0 ) {
 261     daughter = stack->Particle(part->GetFirstDaughter()) ;
 262     decayVx=(Double32_t)daughter->Vx();
 263     decayVy=(Double32_t)daughter->Vy();
 264     decayVz=(Double32_t)daughter->Vz();
 265     decayL = TMath::Sqrt(TMath::Power(partVx-decayVx,2) +
 266                          TMath::Power(partVy-decayVy,2) +
 267                          TMath::Power(partVz-decayVz,2) ) ;
 268     decayRxy = TMath::Sqrt(TMath::Power(decayVx,2) + TMath::Power(decayVy,2) ) ;
 269   }
 270
 271   fCutValues->SetAt(partVx  ,kCutProdVtxXMin) ;
 272   fCutValues->SetAt(partVy  ,kCutProdVtxYMin) ;
 273   fCutValues->SetAt(partVz  ,kCutProdVtxZMin) ;
 274   fCutValues->SetAt(partVx  ,kCutProdVtxXMax) ;
 275   fCutValues->SetAt(partVy  ,kCutProdVtxYMax) ;
 276   fCutValues->SetAt(partVz  ,kCutProdVtxZMax) ;
 277   fCutValues->SetAt(decayVx ,kCutDecVtxXMin)  ;
 278   fCutValues->SetAt(decayVy ,kCutDecVtxYMin)  ;
 279   fCutValues->SetAt(decayVz ,kCutDecVtxZMin)  ;
 280   fCutValues->SetAt(decayVx ,kCutDecVtxXMax)  ;
 281   fCutValues->SetAt(decayVy ,kCutDecVtxYMax)  ;
 282   fCutValues->SetAt(decayVz ,kCutDecVtxZMax)  ;
 283   fCutValues->SetAt(decayL  ,kCutDecLgthMin)  ;
 284   fCutValues->SetAt(decayL  ,kCutDecLgthMax)  ;
 285   fCutValues->SetAt(decayRxy,kCutDecRxyMin)   ;
 286   fCutValues->SetAt(decayRxy,kCutDecRxyMax)   ;
 287
 288   // cut on charge
 289   if ( fRequireIsCharged || fRequireIsNeutral ) {
 290     if (fRequireIsCharged &&  IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
 291     if (fRequireIsNeutral && !IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
 292   }
 293   else fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
 294
 295   // cut on primary/secondary
 296   if ( fRequireIsPrimary || fRequireIsSecondary) {
 297     if (fRequireIsPrimary   &&  IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
 298     if (fRequireIsSecondary && !IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
 299   }
 300   else fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
 301
 302   // cut on PDG code
 303   if ( fRequirePdgCode ) {
 304     if (IsA(mcPart,fPdgCode,kFALSE)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
 305   }
 306   else fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
 307
 308   // cut on decay channel
 309   if ( fDecayChannel ) {
 310     Bool_t goodDecay = kTRUE ;
 311     Short_t nDaughters = mcPart->Particle()->GetNDaughters() ;
 312     if (nDaughters != fDecayChannel->NDaughters()) goodDecay = kFALSE ;
 313     if (goodDecay) {
 314       // now check if decay channel is respected
 315       // first try
 316       for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) {
 317         TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(iDaughter)) ;
 318         if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
 319       }
 320       if (!goodDecay) {
 321         //second try inverting the order of the daughters
 322         goodDecay = kTRUE ;
 323         for (Int_t iDaughter = 0; iDaughter<nDaughters; iDaughter++) {
 324           TParticle* daug = stack->Particle(mcPart->Particle()->GetDaughter(nDaughters-(iDaughter+1))) ;
 325           if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
 326         }
 327       }
 328     }
 329     fCutValues->SetAt((Double32_t)goodDecay,kCutDecayChannel) ;
 330   }
 331   else fCutValues->SetAt((Double32_t)kTRUE,kCutDecayChannel);
 332
 333
 334   // now array of cut is build, fill the bitmap consequently
 335   Int_t iCutBit = -1;
 336   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
 337   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
 338   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
 339   if ( fCutValues->At(++iCutBit) > fProdVtxXMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 340   if ( fCutValues->At(++iCutBit) < fProdVtxXMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 341   if ( fCutValues->At(++iCutBit) > fProdVtxYMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 342   if ( fCutValues->At(++iCutBit) < fProdVtxYMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 343   if ( fCutValues->At(++iCutBit) > fProdVtxZMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 344   if ( fCutValues->At(++iCutBit) < fProdVtxZMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 345   if ( fCutValues->At(++iCutBit) > fDecayVtxXMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 346   if ( fCutValues->At(++iCutBit) < fDecayVtxXMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 347   if ( fCutValues->At(++iCutBit) > fDecayVtxYMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 348   if ( fCutValues->At(++iCutBit) < fDecayVtxYMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 349   if ( fCutValues->At(++iCutBit) > fDecayVtxZMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 350   if ( fCutValues->At(++iCutBit) < fDecayVtxZMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 351   if ( fCutValues->At(++iCutBit) > fDecayLengthMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
 352   if ( fCutValues->At(++iCutBit) < fDecayLengthMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
 353   if ( fCutValues->At(++iCutBit) > fDecayRxyMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 354   if ( fCutValues->At(++iCutBit) < fDecayRxyMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 355   if ( fCutValues->At(++iCutBit) != 0 )             fBitmap->SetBitNumber(iCutBit,kTRUE);
 356 }
 357
 358 //__________________________________________________________________________________
 359 void AliCFParticleGenCuts::SelectionBitMap(AliAODMCParticle* mcPart)
 360 {
 361   //
 362   // test if the track passes the single cuts
 363   // and store the information in a bitmap
 364   //
 365
 366   for (UInt_t i=0; i<kNCuts; i++) {
 367     fBitmap->SetBitNumber(i,kFALSE);
 368     fCutValues->SetAt((Double32_t)0,i) ;
 369   }
 370
 371   // fill the cut array
 372   Double32_t partVx=(Double32_t)mcPart->Xv();
 373   Double32_t partVy=(Double32_t)mcPart->Yv();
 374   Double32_t partVz=(Double32_t)mcPart->Zv();
 375
 376   Double32_t decayVx=0.;
 377   Double32_t decayVy=0.;
 378   Double32_t decayVz=0.;
 379   Double32_t decayL=0.;
 380   Double32_t decayRxy=0.;
 381
 382   AliAODEvent* aod = dynamic_cast<AliAODEvent*>(fMCInfo);
 383   TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
 384   AliAODMCParticle* daughter = 0x0 ;
 385
 386   if (mcPart->GetDaughter(0)>0) {
 387     daughter = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(0)));
 388
 389     decayVx=(Double32_t)daughter->Xv();
 390     decayVy=(Double32_t)daughter->Yv();
 391     decayVz=(Double32_t)daughter->Zv();
 392     decayL = TMath::Sqrt(TMath::Power(partVx-decayVx,2) +
 393                          TMath::Power(partVy-decayVy,2) +
 394                          TMath::Power(partVz-decayVz,2) ) ;
 395     decayRxy = TMath::Sqrt(TMath::Power(decayVx,2) + TMath::Power(decayVy,2) ) ;
 396
 397   }
 398
 399   fCutValues->SetAt(partVx  ,kCutProdVtxXMin) ;
 400   fCutValues->SetAt(partVy  ,kCutProdVtxYMin) ;
 401   fCutValues->SetAt(partVz  ,kCutProdVtxZMin) ;
 402   fCutValues->SetAt(partVx  ,kCutProdVtxXMax) ;
 403   fCutValues->SetAt(partVy  ,kCutProdVtxYMax) ;
 404   fCutValues->SetAt(partVz  ,kCutProdVtxZMax) ;
 405   fCutValues->SetAt(decayVx ,kCutDecVtxXMin)  ;
 406   fCutValues->SetAt(decayVy ,kCutDecVtxYMin)  ;
 407   fCutValues->SetAt(decayVz ,kCutDecVtxZMin)  ;
 408   fCutValues->SetAt(decayVx ,kCutDecVtxXMax)  ;
 409   fCutValues->SetAt(decayVy ,kCutDecVtxYMax)  ;
 410   fCutValues->SetAt(decayVz ,kCutDecVtxZMax)  ;
 411   fCutValues->SetAt(decayL  ,kCutDecLgthMin)  ;
 412   fCutValues->SetAt(decayL  ,kCutDecLgthMax)  ;
 413   fCutValues->SetAt(decayRxy,kCutDecRxyMin)   ;
 414   fCutValues->SetAt(decayRxy,kCutDecRxyMax)   ;
 415
 416   // cut on charge
 417   if ( fRequireIsCharged || fRequireIsNeutral ) {
 418     if (fRequireIsCharged &&  IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
 419     if (fRequireIsNeutral && !IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
 420   }
 421   else fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
 422
 423   // cut on primary/secondary
 424   if ( fRequireIsPrimary || fRequireIsSecondary) {
 425     if (fRequireIsPrimary   &&  IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
 426     if (fRequireIsSecondary && !IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
 427   }
 428   else fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
 429
 430   // cut on PDG code
 431   if ( fRequirePdgCode ) {
 432     if (IsA(mcPart,fPdgCode,kFALSE)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
 433   }
 434   else fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
 435
 436   // cut on decay channel
 437   if ( fDecayChannel ) {
 438     Bool_t goodDecay = kTRUE ;
 439     Short_t nDaughters = 0 ;
 440     if (mcPart->GetDaughter(0) >=0) nDaughters = 1 + mcPart->GetDaughter(1) - mcPart->GetDaughter(0);
 441
 442     if (nDaughters != fDecayChannel->NDaughters()) goodDecay = kFALSE ;
 443     if (goodDecay) {
 444       // now check if decay channel is respected
 445       // first try
 446       for (Int_t iDaughter = 0 ; iDaughter<nDaughters; iDaughter++) {
 447         AliAODMCParticle* daug = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(0)+iDaughter));
 448         if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
 449       }
 450       if (!goodDecay) {
 451         //second try inverting the order of the daughters
 452         goodDecay = kTRUE ;
 453         for (Int_t iDaughter = 0 ; iDaughter<nDaughters; iDaughter++) {
 454           AliAODMCParticle* daug = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(1)-iDaughter));
 455           if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
 456         }
 457       }
 458     }
 459     fCutValues->SetAt((Double32_t)goodDecay,kCutDecayChannel) ;
 460   }
 461   else fCutValues->SetAt((Double32_t)kTRUE,kCutDecayChannel);
 462
 463
 464   // now array of cut is build, fill the bitmap consequently
 465   Int_t iCutBit = -1;
 466   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
 467   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
 468   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
 469   if ( fCutValues->At(++iCutBit) > fProdVtxXMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 470   if ( fCutValues->At(++iCutBit) < fProdVtxXMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 471   if ( fCutValues->At(++iCutBit) > fProdVtxYMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 472   if ( fCutValues->At(++iCutBit) < fProdVtxYMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 473   if ( fCutValues->At(++iCutBit) > fProdVtxZMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 474   if ( fCutValues->At(++iCutBit) < fProdVtxZMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 475   if ( fCutValues->At(++iCutBit) > fDecayVtxXMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 476   if ( fCutValues->At(++iCutBit) < fDecayVtxXMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 477   if ( fCutValues->At(++iCutBit) > fDecayVtxYMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 478   if ( fCutValues->At(++iCutBit) < fDecayVtxYMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 479   if ( fCutValues->At(++iCutBit) > fDecayVtxZMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 480   if ( fCutValues->At(++iCutBit) < fDecayVtxZMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 481   if ( fCutValues->At(++iCutBit) > fDecayLengthMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
 482   if ( fCutValues->At(++iCutBit) < fDecayLengthMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
 483   if ( fCutValues->At(++iCutBit) > fDecayRxyMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 484   if ( fCutValues->At(++iCutBit) < fDecayRxyMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 485   if ( fCutValues->At(++iCutBit) != 0 )             fBitmap->SetBitNumber(iCutBit,kTRUE);
 486 }
 487
 488
 489 //__________________________________________________________________________________
 490 void AliCFParticleGenCuts::FillHistograms(TObject* /*obj*/, Bool_t afterCuts)
 491 {
 492   //
 493   // fill the QA histograms
 494   //
 495
 496   for (int iCutNumber = 0; iCutNumber < kNCuts; iCutNumber++)
 497     fhQA[iCutNumber][afterCuts]->Fill(fCutValues->At(iCutNumber));
 498
 499   // fill cut statistics and cut correlation histograms with information from the bitmap
 500   if (afterCuts) return;
 501
 502   // Number of single cuts in this class
 503   UInt_t ncuts = fBitmap->GetNbits();
 504   for(UInt_t bit=0; bit<ncuts;bit++) {
 505     if (!fBitmap->TestBitNumber(bit)) {
 506       fhCutStatistics->Fill(bit+1);
 507       for (UInt_t bit2=bit; bit2<ncuts;bit2++) {
 508         if (!fBitmap->TestBitNumber(bit2))
 509           fhCutCorrelation->Fill(bit+1,bit2+1);
 510       }
 511     }
 512   }
 513 }
 514
 515 //__________________________________________________________________________________
 516 void AliCFParticleGenCuts::AddQAHistograms(TList *qaList) {
 517   //
 518   // saves the histograms in a TList
 519   //
 520
 521   DefineHistograms();
 522
 523   qaList->Add(fhCutStatistics);
 524   qaList->Add(fhCutCorrelation);
 525
 526   for (Int_t j=0; j<kNStepQA; j++) {
 527     for(Int_t i=0; i<kNCuts; i++)
 528       qaList->Add(fhQA[i][j]);
 529   }
 530 }
 531
 532 //__________________________________________________________________________________
 533 void AliCFParticleGenCuts::DefineHistograms() {
 534   //
 535   // histograms for cut variables, cut statistics and cut correlations
 536   //
 537   Int_t color = 2;
 538
 539   // book cut statistics and cut correlation histograms
 540   fhCutStatistics = new TH1F(Form("%s_cut_statistics",GetName()),"",kNCuts,0.5,kNCuts+0.5);
 541   fhCutStatistics->SetLineWidth(2);
 542   int k = 1;
 543   fhCutStatistics->GetXaxis()->SetBinLabel(k,"charge")     ; k++;
 544   fhCutStatistics->GetXaxis()->SetBinLabel(k,"prim/sec")   ; k++;
 545   fhCutStatistics->GetXaxis()->SetBinLabel(k,"PDG")        ; k++;
 546   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxXMin")    ; k++;
 547   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxXMax")    ; k++;
 548   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxYMin")    ; k++;
 549   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxYMax")    ; k++;
 550   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxZMin")    ; k++;
 551   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMax")    ; k++;
 552   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecXMin")    ; k++;
 553   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecXMax")    ; k++;
 554   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecYMin")    ; k++;
 555   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecYMax")    ; k++;
 556   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMin")    ; k++;
 557   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMax")    ; k++;
 558   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecLgthMin") ; k++;
 559   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecLgthMax") ; k++;
 560   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMin")  ; k++;
 561   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMax")  ; k++;
 562   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecChannel") ; k++;
 563
 564
 565   fhCutCorrelation = new TH2F(Form("%s_cut_correlation",GetName()),"",kNCuts,0.5,kNCuts+0.5,kNCuts,0.5,kNCuts+0.5);
 566   fhCutCorrelation->SetLineWidth(2);
 567   for (k=1; k<=kNCuts; k++) {
 568     fhCutCorrelation->GetXaxis()->SetBinLabel(k,fhCutStatistics->GetXaxis()->GetBinLabel(k));
 569     fhCutCorrelation->GetYaxis()->SetBinLabel(k,fhCutStatistics->GetXaxis()->GetBinLabel(k));
 570   }
 571
 572   Char_t str[256];
 573   for (int i=0; i<kNStepQA; i++) {
 574     if (i==0) sprintf(str," ");
 575     else sprintf(str,"_cut");
 576     fhQA[kCutCharge]      [i] = new TH1F(Form("%s_charge%s"      ,GetName(),str),"",2,0,2);
 577     fhQA[kCutPrimSec]     [i] = new TH1F(Form("%s_primSec%s"     ,GetName(),str),"",2,0,2);
 578     fhQA[kCutPDGCode]     [i] = new TH1F(Form("%s_pdgCode%s"     ,GetName(),str),"",2,0,2);
 579     fhQA[kCutProdVtxXMin] [i] = new TH1F(Form("%s_prodVtxXMin%s" ,GetName(),str),"",100,0,10);
 580     fhQA[kCutProdVtxXMax] [i] = new TH1F(Form("%s_prodVtxXMax%s" ,GetName(),str),"",100,0,10);
 581     fhQA[kCutProdVtxYMin] [i] = new TH1F(Form("%s_prodVtxYMin%s" ,GetName(),str),"",100,0,10);
 582     fhQA[kCutProdVtxYMax] [i] = new TH1F(Form("%s_prodVtxYMax%s" ,GetName(),str),"",100,0,10);
 583     fhQA[kCutProdVtxZMin] [i] = new TH1F(Form("%s_prodVtxZMin%s" ,GetName(),str),"",100,0,10);
 584     fhQA[kCutProdVtxZMax] [i] = new TH1F(Form("%s_prodVtxZMax%s" ,GetName(),str),"",100,0,10);
 585     fhQA[kCutDecVtxXMin]  [i] = new TH1F(Form("%s_decVtxXMin%s"  ,GetName(),str),"",100,0,10);
 586     fhQA[kCutDecVtxXMax]  [i] = new TH1F(Form("%s_decVtxXMax%s"  ,GetName(),str),"",100,0,10);
 587     fhQA[kCutDecVtxYMin]  [i] = new TH1F(Form("%s_decVtxYMin%s"  ,GetName(),str),"",100,0,10);
 588     fhQA[kCutDecVtxYMax]  [i] = new TH1F(Form("%s_decVtxYMax%s"  ,GetName(),str),"",100,0,10);
 589     fhQA[kCutDecVtxZMin]  [i] = new TH1F(Form("%s_decVtxZMin%s"  ,GetName(),str),"",100,0,10);
 590     fhQA[kCutDecVtxZMax]  [i] = new TH1F(Form("%s_decVtxZMax%s"  ,GetName(),str),"",100,0,10);
 591     fhQA[kCutDecLgthMin]  [i] = new TH1F(Form("%s_decLengthMin%s",GetName(),str),"",100,0,10);
 592     fhQA[kCutDecLgthMax]  [i] = new TH1F(Form("%s_decLengthMax%s",GetName(),str),"",100,0,10);
 593     fhQA[kCutDecRxyMin]   [i] = new TH1F(Form("%s_decRxyMin%s"   ,GetName(),str),"",100,0,10);
 594     fhQA[kCutDecRxyMax]   [i] = new TH1F(Form("%s_decRxyMax%s"   ,GetName(),str),"",100,0,10);
 595     fhQA[kCutDecayChannel][i] = new TH1F(Form("%s_decayChannel%s",GetName(),str),"",2,0,2);
 596   }
 597   for(Int_t i=0; i<kNCuts; i++) fhQA[i][1]->SetLineColor(color);
 598 }
 599
 600
 601 //______________________________
 602 Bool_t AliCFParticleGenCuts::IsCharged(AliVParticle *mcPart) {
 603   //
 604   //check if particle is charged.
 605   //
 606   if (mcPart->Charge()==0) return kFALSE;
 607   return kTRUE;
 608 }
 609 //______________________________
 610 Bool_t AliCFParticleGenCuts::IsPrimary(AliMCParticle *mcPart) {
 611   //
 612   //check if particle is primary (standard definition)
 613   //
 614
 615   AliStack* stack = ((AliMCEvent*)fMCInfo)->Stack();
 616
 617   if (!stack->IsPhysicalPrimary(mcPart->GetLabel())) return kFALSE;
 618   return kTRUE;
 619 }
 620 //______________________________
 621 Bool_t AliCFParticleGenCuts::IsPrimary(AliAODMCParticle *mcPart) {
 622   //
 623   //check if particle is primary (standard definition)
 624   //
 625
 626   if (!mcPart->IsPhysicalPrimary()) return kFALSE;
 627   return kTRUE;
 628 }
 629 //______________________________
 630 Bool_t AliCFParticleGenCuts::IsPrimaryCharged(AliVParticle *mcPart) {
 631   //
 632   //check if a charged particle is primary (standard definition)
 633   //
 634
 635   if (!fIsAODMC) {
 636     if (!IsPrimary((AliMCParticle*)mcPart) || !IsCharged(mcPart)) return kFALSE ;
 637   }
 638   else {
 639     if (!IsPrimary((AliAODMCParticle*)mcPart) || !IsCharged(mcPart)) return kFALSE ;
 640   }
 641   return kTRUE;
 642 }
 643 //______________________________
 644 Bool_t AliCFParticleGenCuts::IsA(AliMCParticle *mcPart, Int_t pdg, Bool_t abs) {
 645   //
 646   //Check on the pdg code of the MC particle. if abs=kTRUE then check on the
 647   //absolute value. By default is set to kFALSE.
 648   //
 649   TParticle* part = mcPart->Particle();
 650   Int_t pdgCode = part->GetPdgCode();
 651
 652   if (abs) {
 653     pdgCode = TMath::Abs(pdgCode);
 654     pdg = TMath::Abs(pdg);
 655   }
 656   if (pdgCode != pdg ) return kFALSE;
 657   return kTRUE;
 658 }
 659 //______________________________
 660 Bool_t AliCFParticleGenCuts::IsA(AliAODMCParticle *mcPart, Int_t pdg, Bool_t abs) {
 661   //
 662   //Check on the pdg code of the MC particle. if abs=kTRUE then check on the
 663   //absolute value. By default is set to kFALSE.
 664   //
 665   Int_t pdgCode = mcPart->GetPdgCode();
 666
 667   if (abs) {
 668     pdgCode = TMath::Abs(pdgCode);
 669     pdg = TMath::Abs(pdg);
 670   }
 671   if (pdgCode != pdg ) return kFALSE;
 672   return kTRUE;
 673 }
 674 //______________________________
 675 void AliCFParticleGenCuts::SetEvtInfo(TObject* mcEvent) {
 676   //
 677   // Sets pointer to MC event information (AliMCEvent)
 678   //
 679
 680   if (!mcEvent) {
 681     AliError("Pointer to MC Event is null !");
 682     return;
 683   }
 684
 685   TString className(mcEvent->ClassName());
 686   if (className.CompareTo("AliMCEvent") != 0 && className.CompareTo("AliAODEvent") != 0) {
 687     AliError("argument must point to an AliMCEvent or an AliAODEvent !");
 688     return ;
 689   }
 690
 691   if (fIsAODMC) fMCInfo = dynamic_cast<AliAODEvent*>(mcEvent) ;
 692   else          fMCInfo = dynamic_cast<AliMCEvent*> (mcEvent) ;
 693 }