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