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 && prodVtxXYmin >= 1))     fBitmap->SetBitNumber(iCutBit,kTRUE);
 389
 390   ++iCutBit;
 391   if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxXMax)
 392     || ( fProdVtxRange2D && prodVtxXYmax <= 1))     fBitmap->SetBitNumber(iCutBit,kTRUE);
 393
 394   ++iCutBit;
 395   if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxYMin)
 396     || ( fProdVtxRange2D && prodVtxXYmin >= 1))     fBitmap->SetBitNumber(iCutBit,kTRUE);
 397
 398   ++iCutBit;
 399   if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxYMax)
 400     || ( fProdVtxRange2D && prodVtxXYmax <= 1))     fBitmap->SetBitNumber(iCutBit,kTRUE);
 401
 402   if ( fCutValues->At(++iCutBit) > fProdVtxZMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 403   if ( fCutValues->At(++iCutBit) < fProdVtxZMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 404   if ( fCutValues->At(++iCutBit) > fDecayVtxXMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 405   if ( fCutValues->At(++iCutBit) < fDecayVtxXMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 406   if ( fCutValues->At(++iCutBit) > fDecayVtxYMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 407   if ( fCutValues->At(++iCutBit) < fDecayVtxYMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 408   if ( fCutValues->At(++iCutBit) > fDecayVtxZMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 409   if ( fCutValues->At(++iCutBit) < fDecayVtxZMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 410   if ( fCutValues->At(++iCutBit) > fDecayLengthMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
 411   if ( fCutValues->At(++iCutBit) < fDecayLengthMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
 412   if ( fCutValues->At(++iCutBit) > fDecayRxyMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 413   if ( fCutValues->At(++iCutBit) < fDecayRxyMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 414   if ( fCutValues->At(++iCutBit) != 0 )             fBitmap->SetBitNumber(iCutBit,kTRUE);
 415 }
 416
 417 //__________________________________________________________________________________
 418 void AliCFParticleGenCuts::SelectionBitMap(AliAODMCParticle* mcPart)
 419 {
 420   //
 421   // test if the track passes the single cuts
 422   // and store the information in a bitmap
 423   //
 424
 425   for (UInt_t i=0; i<kNCuts; i++) {
 426     fBitmap->SetBitNumber(i,kFALSE);
 427     fCutValues->SetAt((Double32_t)0,i) ;
 428   }
 429
 430   // fill the cut array
 431   Double32_t partVx=(Double32_t)mcPart->Xv();
 432   Double32_t partVy=(Double32_t)mcPart->Yv();
 433   Double32_t partVz=(Double32_t)mcPart->Zv();
 434
 435   // calculate the production vertex ellipse
 436   Double32_t prodVtxXYmin = 0.;
 437   if (fProdVtxXMin!=0 && fProdVtxYMin!=0)
 438         prodVtxXYmin = partVx*partVx/(fProdVtxXMin*fProdVtxXMin) + partVy*partVy/(fProdVtxYMin*fProdVtxYMin);
 439   Double32_t prodVtxXYmax = 0.;
 440   if(fProdVtxXMax!=0 && fProdVtxYMax!=0)
 441         prodVtxXYmax = partVx*partVx/(fProdVtxXMax*fProdVtxXMax) + partVy*partVy/(fProdVtxYMax*fProdVtxYMax);
 442
 443   Double32_t decayVx=0.;
 444   Double32_t decayVy=0.;
 445   Double32_t decayVz=0.;
 446   Double32_t decayL=0.;
 447   Double32_t decayRxy=0.;
 448
 449   AliAODEvent* aod = dynamic_cast<AliAODEvent*>(fMCInfo);
 450   TClonesArray* mcArray = dynamic_cast<TClonesArray*>(aod->FindListObject(AliAODMCParticle::StdBranchName()));
 451   AliAODMCParticle* daughter = 0x0 ;
 452
 453   if (mcPart->GetDaughter(0)>0) {
 454     daughter = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(0)));
 455
 456     decayVx=(Double32_t)daughter->Xv();
 457     decayVy=(Double32_t)daughter->Yv();
 458     decayVz=(Double32_t)daughter->Zv();
 459     decayL = TMath::Sqrt(TMath::Power(partVx-decayVx,2) +
 460                          TMath::Power(partVy-decayVy,2) +
 461                          TMath::Power(partVz-decayVz,2) ) ;
 462     decayRxy = TMath::Sqrt(TMath::Power(decayVx,2) + TMath::Power(decayVy,2) ) ;
 463
 464   }
 465
 466   fCutValues->SetAt(partVx  ,kCutProdVtxXMin) ;
 467   fCutValues->SetAt(partVy  ,kCutProdVtxYMin) ;
 468   fCutValues->SetAt(partVz  ,kCutProdVtxZMin) ;
 469   fCutValues->SetAt(partVx  ,kCutProdVtxXMax) ;
 470   fCutValues->SetAt(partVy  ,kCutProdVtxYMax) ;
 471   fCutValues->SetAt(partVz  ,kCutProdVtxZMax) ;
 472   fCutValues->SetAt(decayVx ,kCutDecVtxXMin)  ;
 473   fCutValues->SetAt(decayVy ,kCutDecVtxYMin)  ;
 474   fCutValues->SetAt(decayVz ,kCutDecVtxZMin)  ;
 475   fCutValues->SetAt(decayVx ,kCutDecVtxXMax)  ;
 476   fCutValues->SetAt(decayVy ,kCutDecVtxYMax)  ;
 477   fCutValues->SetAt(decayVz ,kCutDecVtxZMax)  ;
 478   fCutValues->SetAt(decayL  ,kCutDecLgthMin)  ;
 479   fCutValues->SetAt(decayL  ,kCutDecLgthMax)  ;
 480   fCutValues->SetAt(decayRxy,kCutDecRxyMin)   ;
 481   fCutValues->SetAt(decayRxy,kCutDecRxyMax)   ;
 482
 483   // cut on charge
 484   if ( fRequireIsCharged || fRequireIsNeutral ) {
 485     if (fRequireIsCharged &&  IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
 486     if (fRequireIsNeutral && !IsCharged(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
 487   }
 488   else fCutValues->SetAt((Double32_t)kTRUE,kCutCharge) ;
 489
 490   // cut on primary/secondary
 491   if ( fRequireIsPrimary || fRequireIsSecondary) {
 492     if (fRequireIsPrimary   &&  IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
 493     if (fRequireIsSecondary && !IsPrimary(mcPart)) fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
 494   }
 495   else fCutValues->SetAt((Double32_t)kTRUE,kCutPrimSec);
 496
 497   // cut on PDG code
 498   if ( fRequirePdgCode ) {
 499     if (IsA(mcPart,fPdgCode,fRequireAbsolutePdg)) fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
 500   }
 501   else fCutValues->SetAt((Double32_t)kTRUE,kCutPDGCode);
 502
 503   // cut on decay channel
 504   if ( fDecayChannel ) {
 505     Bool_t goodDecay = kTRUE ;
 506     Short_t nDaughters = 0 ;
 507     if (mcPart->GetDaughter(0) >=0) nDaughters = 1 + mcPart->GetDaughter(1) - mcPart->GetDaughter(0);
 508
 509     if (nDaughters != fDecayChannel->NDaughters()) goodDecay = kFALSE ;
 510     if (goodDecay) {
 511       // now check if decay channel is respected
 512       // first try
 513       for (Int_t iDaughter = 0 ; iDaughter<nDaughters; iDaughter++) {
 514         AliAODMCParticle* daug = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(0)+iDaughter));
 515         if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
 516       }
 517       if (!goodDecay) {
 518         //second try inverting the order of the daughters
 519         goodDecay = kTRUE ;
 520         for (Int_t iDaughter = 0 ; iDaughter<nDaughters; iDaughter++) {
 521           AliAODMCParticle* daug = dynamic_cast<AliAODMCParticle*>(mcArray->At(mcPart->GetDaughter(1)-iDaughter));
 522           if (daug->GetPdgCode() != fDecayChannel->DaughterPdgCode(iDaughter)) {goodDecay = kFALSE; break;}
 523         }
 524       }
 525     }
 526     fCutValues->SetAt((Double32_t)goodDecay,kCutDecayChannel) ;
 527   }
 528   else fCutValues->SetAt((Double32_t)kTRUE,kCutDecayChannel);
 529
 530
 531   // now array of cut is build, fill the bitmap consequently
 532   Int_t iCutBit = -1;
 533   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
 534   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
 535   if ( fCutValues->At(++iCutBit) !=0 )              fBitmap->SetBitNumber(iCutBit,kTRUE);
 536
 537   ++iCutBit;
 538   if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxXMin)
 539     || ( fProdVtxRange2D && prodVtxXYmin >= 1))     fBitmap->SetBitNumber(iCutBit,kTRUE);
 540
 541   ++iCutBit;
 542   if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxXMax)
 543     || ( fProdVtxRange2D && prodVtxXYmax <= 1))     fBitmap->SetBitNumber(iCutBit,kTRUE);
 544
 545   ++iCutBit;
 546   if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) > fProdVtxYMin)
 547     || ( fProdVtxRange2D && prodVtxXYmin >= 1))     fBitmap->SetBitNumber(iCutBit,kTRUE);
 548
 549   ++iCutBit;
 550   if ( (!fProdVtxRange2D && fCutValues->At(iCutBit) < fProdVtxYMax)
 551     || ( fProdVtxRange2D && prodVtxXYmax <= 1))     fBitmap->SetBitNumber(iCutBit,kTRUE);
 552
 553   if ( fCutValues->At(++iCutBit) > fProdVtxZMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 554   if ( fCutValues->At(++iCutBit) < fProdVtxZMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 555   if ( fCutValues->At(++iCutBit) > fDecayVtxXMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 556   if ( fCutValues->At(++iCutBit) < fDecayVtxXMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 557   if ( fCutValues->At(++iCutBit) > fDecayVtxYMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 558   if ( fCutValues->At(++iCutBit) < fDecayVtxYMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 559   if ( fCutValues->At(++iCutBit) > fDecayVtxZMin)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 560   if ( fCutValues->At(++iCutBit) < fDecayVtxZMax)   fBitmap->SetBitNumber(iCutBit,kTRUE);
 561   if ( fCutValues->At(++iCutBit) > fDecayLengthMin) fBitmap->SetBitNumber(iCutBit,kTRUE);
 562   if ( fCutValues->At(++iCutBit) < fDecayLengthMax) fBitmap->SetBitNumber(iCutBit,kTRUE);
 563   if ( fCutValues->At(++iCutBit) > fDecayRxyMin)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 564   if ( fCutValues->At(++iCutBit) < fDecayRxyMax)    fBitmap->SetBitNumber(iCutBit,kTRUE);
 565   if ( fCutValues->At(++iCutBit) != 0 )             fBitmap->SetBitNumber(iCutBit,kTRUE);
 566 }
 567
 568
 569 //__________________________________________________________________________________
 570 void AliCFParticleGenCuts::FillHistograms(TObject* /*obj*/, Bool_t afterCuts)
 571 {
 572   //
 573   // fill the QA histograms
 574   //
 575
 576   for (int iCutNumber = 0; iCutNumber < kNCuts; iCutNumber++)
 577     fhQA[iCutNumber][afterCuts]->Fill(fCutValues->At(iCutNumber));
 578
 579   fhProdVtxXY[afterCuts]->Fill(fCutValues->At(4),fCutValues->At(5));
 580
 581   // fill cut statistics and cut correlation histograms with information from the bitmap
 582   if (afterCuts) return;
 583
 584   // Number of single cuts in this class
 585   UInt_t ncuts = fBitmap->GetNbits();
 586   for(UInt_t bit=0; bit<ncuts;bit++) {
 587     if (!fBitmap->TestBitNumber(bit)) {
 588       fhCutStatistics->Fill(bit+1);
 589       for (UInt_t bit2=bit; bit2<ncuts;bit2++) {
 590         if (!fBitmap->TestBitNumber(bit2))
 591           fhCutCorrelation->Fill(bit+1,bit2+1);
 592       }
 593     }
 594   }
 595 }
 596
 597 //__________________________________________________________________________________
 598 void AliCFParticleGenCuts::AddQAHistograms(TList *qaList) {
 599   //
 600   // saves the histograms in a TList
 601   //
 602
 603   DefineHistograms();
 604
 605   qaList->Add(fhCutStatistics);
 606   qaList->Add(fhCutCorrelation);
 607
 608   for (Int_t j=0; j<kNStepQA; j++) {
 609     qaList->Add(fhProdVtxXY[j]);
 610     for(Int_t i=0; i<kNCuts; i++)
 611       qaList->Add(fhQA[i][j]);
 612   }
 613 }
 614
 615 //__________________________________________________________________________________
 616 void AliCFParticleGenCuts::DefineHistograms() {
 617   //
 618   // histograms for cut variables, cut statistics and cut correlations
 619   //
 620   Int_t color = 2;
 621
 622   // book cut statistics and cut correlation histograms
 623   fhCutStatistics = new TH1F(Form("%s_cut_statistics",GetName()),"",kNCuts,0.5,kNCuts+0.5);
 624   fhCutStatistics->SetLineWidth(2);
 625   int k = 1;
 626   fhCutStatistics->GetXaxis()->SetBinLabel(k,"charge")     ; k++;
 627   fhCutStatistics->GetXaxis()->SetBinLabel(k,"prim/sec")   ; k++;
 628   fhCutStatistics->GetXaxis()->SetBinLabel(k,"PDG")        ; k++;
 629   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxXMin")    ; k++;
 630   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxXMax")    ; k++;
 631   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxYMin")    ; k++;
 632   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxYMax")    ; k++;
 633   fhCutStatistics->GetXaxis()->SetBinLabel(k,"VtxZMin")    ; k++;
 634   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMax")    ; k++;
 635   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecXMin")    ; k++;
 636   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecXMax")    ; k++;
 637   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecYMin")    ; k++;
 638   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecYMax")    ; k++;
 639   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMin")    ; k++;
 640   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecZMax")    ; k++;
 641   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecLgthMin") ; k++;
 642   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecLgthMax") ; k++;
 643   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMin")  ; k++;
 644   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecRxyMax")  ; k++;
 645   fhCutStatistics->GetXaxis()->SetBinLabel(k,"DecChannel") ; k++;
 646
 647
 648   fhCutCorrelation = new TH2F(Form("%s_cut_correlation",GetName()),"",kNCuts,0.5,kNCuts+0.5,kNCuts,0.5,kNCuts+0.5);
 649   fhCutCorrelation->SetLineWidth(2);
 650   for (k=1; k<=kNCuts; k++) {
 651     fhCutCorrelation->GetXaxis()->SetBinLabel(k,fhCutStatistics->GetXaxis()->GetBinLabel(k));
 652     fhCutCorrelation->GetYaxis()->SetBinLabel(k,fhCutStatistics->GetXaxis()->GetBinLabel(k));
 653   }
 654
 655   Char_t str[256];
 656   for (int i=0; i<kNStepQA; i++) {
 657     if (i==0) sprintf(str," ");
 658     else sprintf(str,"_cut");
 659     fhQA[kCutCharge]      [i] = new TH1F(Form("%s_charge%s"      ,GetName(),str),"",2,0,2);
 660     fhQA[kCutPrimSec]     [i] = new TH1F(Form("%s_primSec%s"     ,GetName(),str),"",2,0,2);
 661     fhQA[kCutPDGCode]     [i] = new TH1F(Form("%s_pdgCode%s"     ,GetName(),str),"",2,0,2);
 662     fhQA[kCutProdVtxXMin] [i] = new TH1F(Form("%s_prodVtxXMin%s" ,GetName(),str),"",100,-10,10);
 663     fhQA[kCutProdVtxXMax] [i] = new TH1F(Form("%s_prodVtxXMax%s" ,GetName(),str),"",100,-10,10);
 664     fhQA[kCutProdVtxYMin] [i] = new TH1F(Form("%s_prodVtxYMin%s" ,GetName(),str),"",100,-10,10);
 665     fhQA[kCutProdVtxYMax] [i] = new TH1F(Form("%s_prodVtxYMax%s" ,GetName(),str),"",100,-10,10);
 666     fhQA[kCutProdVtxZMin] [i] = new TH1F(Form("%s_prodVtxZMin%s" ,GetName(),str),"",100,-10,10);
 667     fhQA[kCutProdVtxZMax] [i] = new TH1F(Form("%s_prodVtxZMax%s" ,GetName(),str),"",100,-10,10);
 668     fhQA[kCutDecVtxXMin]  [i] = new TH1F(Form("%s_decVtxXMin%s"  ,GetName(),str),"",100,0,10);
 669     fhQA[kCutDecVtxXMax]  [i] = new TH1F(Form("%s_decVtxXMax%s"  ,GetName(),str),"",100,0,10);
 670     fhQA[kCutDecVtxYMin]  [i] = new TH1F(Form("%s_decVtxYMin%s"  ,GetName(),str),"",100,0,10);
 671     fhQA[kCutDecVtxYMax]  [i] = new TH1F(Form("%s_decVtxYMax%s"  ,GetName(),str),"",100,0,10);
 672     fhQA[kCutDecVtxZMin]  [i] = new TH1F(Form("%s_decVtxZMin%s"  ,GetName(),str),"",100,0,10);
 673     fhQA[kCutDecVtxZMax]  [i] = new TH1F(Form("%s_decVtxZMax%s"  ,GetName(),str),"",100,0,10);
 674     fhQA[kCutDecLgthMin]  [i] = new TH1F(Form("%s_decLengthMin%s",GetName(),str),"",100,0,10);
 675     fhQA[kCutDecLgthMax]  [i] = new TH1F(Form("%s_decLengthMax%s",GetName(),str),"",100,0,10);
 676     fhQA[kCutDecRxyMin]   [i] = new TH1F(Form("%s_decRxyMin%s"   ,GetName(),str),"",100,0,10);
 677     fhQA[kCutDecRxyMax]   [i] = new TH1F(Form("%s_decRxyMax%s"   ,GetName(),str),"",100,0,10);
 678     fhQA[kCutDecayChannel][i] = new TH1F(Form("%s_decayChannel%s",GetName(),str),"",2,0,2);
 679     fhProdVtxXY           [i] = new TH2F(Form("%s_prodVtxXY%s"   ,GetName(),str),"",100,0,10,100,0,10);
 680     fhProdVtxXY           [i] ->GetXaxis()->SetTitle("x_{production vertex}");
 681     fhProdVtxXY           [i] ->GetYaxis()->SetTitle("y_{production vertex}");
 682     fhQA[kCutProdVtxXMax] [i] ->GetXaxis()->SetTitle("x_{production vertex}");
 683     fhQA[kCutProdVtxYMax] [i] ->GetXaxis()->SetTitle("y_{production vertex}");
 684   }
 685   for(Int_t i=0; i<kNCuts; i++) fhQA[i][1]->SetLineColor(color);
 686 }
 687
 688
 689 //______________________________
 690 Bool_t AliCFParticleGenCuts::IsCharged(AliVParticle *mcPart) {
 691   //
 692   //check if particle is charged.
 693   //
 694   if (mcPart->Charge()==0) return kFALSE;
 695   return kTRUE;
 696 }
 697 //______________________________
 698 Bool_t AliCFParticleGenCuts::IsPrimary(AliMCParticle *mcPart) {
 699   //
 700   //check if particle is primary (standard definition)
 701   //
 702
 703   AliStack* stack = ((AliMCEvent*)fMCInfo)->Stack();
 704
 705   if (!stack->IsPhysicalPrimary(mcPart->GetLabel())) return kFALSE;
 706   return kTRUE;
 707 }
 708 //______________________________
 709 Bool_t AliCFParticleGenCuts::IsPrimary(AliAODMCParticle *mcPart) {
 710   //
 711   //check if particle is primary (standard definition)
 712   //
 713
 714   if (!mcPart->IsPhysicalPrimary()) return kFALSE;
 715   return kTRUE;
 716 }
 717 //______________________________
 718 Bool_t AliCFParticleGenCuts::IsPrimaryCharged(AliVParticle *mcPart) {
 719   //
 720   //check if a charged particle is primary (standard definition)
 721   //
 722
 723   if (!fIsAODMC) {
 724     if (!IsPrimary((AliMCParticle*)mcPart) || !IsCharged(mcPart)) return kFALSE ;
 725   }
 726   else {
 727     if (!IsPrimary((AliAODMCParticle*)mcPart) || !IsCharged(mcPart)) return kFALSE ;
 728   }
 729   return kTRUE;
 730 }
 731 //______________________________
 732 Bool_t AliCFParticleGenCuts::IsA(AliMCParticle *mcPart, Int_t pdg, Bool_t abs) {
 733   //
 734   //Check on the pdg code of the MC particle. if abs=kTRUE then check on the
 735   //absolute value.
 736   //
 737   TParticle* part = mcPart->Particle();
 738   Int_t pdgCode = part->GetPdgCode();
 739
 740   if (abs) {
 741     pdgCode = TMath::Abs(pdgCode);
 742     pdg = TMath::Abs(pdg);
 743   }
 744   if (pdgCode != pdg ) return kFALSE;
 745   return kTRUE;
 746 }
 747 //______________________________
 748 Bool_t AliCFParticleGenCuts::IsA(AliAODMCParticle *mcPart, Int_t pdg, Bool_t abs) {
 749   //
 750   //Check on the pdg code of the MC particle. if abs=kTRUE then check on the
 751   //absolute value.
 752   //
 753   Int_t pdgCode = mcPart->GetPdgCode();
 754
 755   if (abs) {
 756     pdgCode = TMath::Abs(pdgCode);
 757     pdg = TMath::Abs(pdg);
 758   }
 759   if (pdgCode != pdg ) return kFALSE;
 760   return kTRUE;
 761 }
 762 //______________________________
 763 void AliCFParticleGenCuts::SetMCEventInfo(const TObject* mcEvent) {
 764   //
 765   // Sets pointer to MC event information (AliMCEvent)
 766   //
 767
 768   if (!mcEvent) {
 769     AliError("Pointer to MC Event is null !");
 770     return;
 771   }
 772
 773   TString className(mcEvent->ClassName());
 774   if (className.CompareTo("AliMCEvent") != 0 && className.CompareTo("AliAODEvent") != 0) {
 775     AliError("argument must point to an AliMCEvent or an AliAODEvent !");
 776     return ;
 777   }
 778
 779   fMCInfo = (AliVEvent*)mcEvent ;
 780 }