PWGDQ/dielectron/AliDielectronV0Cuts.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2009, 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 //   Cut class providing cuts to V0 candidates                           //
  18 //   Selection or deselection of V0 candiates can be done.               //
  19 //                                                                       //
  20 // Authors:                                                              //
  21 //   Julian Book <Julian.Book@cern.ch>                                   //
  22 /*
  23
  24 Class to provide some V0 selection using exactley the same cuts in ESDs as in AODs.
  25 It implements the PID cut class AliDielectronPID and the standard AliDielectronVarCuts for
  26 the configuration of leg respective pair cuts. These pair cuts are applied on the KFparticle
  27 build by the legs.
  28
  29 Some QA cuts for the tracks are applied before the V0 pair is build. These cuts are:
  30 AliDielectronVarCuts dauQAcuts1;
  31 dauQAcuts1.AddCut(AliDielectronVarManager::kPt,            0.05, 100.0);
  32 dauQAcuts1.AddCut(AliDielectronVarManager::kEta,          -0.9,    0.9);
  33 dauQAcuts1.AddCut(AliDielectronVarManager::kNclsTPC,      50.0,  160.0);
  34 dauQAcuts1.AddCut(AliDielectronVarManager::kTPCchi2Cl,     0.0,    4.0);
  35 AliDielectronTrackCuts dauQAcuts2;
  36 dauQAcuts2.SetRequireTPCRefit(kTRUE);
  37
  38
  39
  40 Example configuration:
  41
  42   AliDielectronV0Cuts *gammaV0Cuts = new AliDielectronV0Cuts("IsGamma","IsGamma");
  43
  44   // which V0 finder you want to use
  45   gammaV0Cuts->SetV0finder(kOnTheFly);  // kAll(default), kOffline or kOnTheFly
  46
  47   // add some pdg codes (they are used then by the KF package and important for gamma conversions)
  48   gammaV0Cuts->SetPdgCodes(22,11,11); // mother, daughter1 and 2
  49
  50   // add default PID cuts (defined in AliDielectronPID)
  51   // requirement can be set to at least one(kAny) of the tracks or to both(kBoth)
  52   gammaV0Cuts->SetDefaultPID(16, AliDielectronV0Cuts::kAny);
  53
  54   // add the pair cuts for V0 candidates
  55   gammaV0Cuts->AddCut(AliDielectronVarManager::kCosPointingAngle, TMath::Cos(0.02),   1.0, kFALSE);
  56   gammaV0Cuts->AddCut(AliDielectronVarManager::kChi2NDF,                       0.0,  10.0, kFALSE);
  57   gammaV0Cuts->AddCut(AliDielectronVarManager::kLegDist,                       0.0,   0.25, kFALSE);
  58   gammaV0Cuts->AddCut(AliDielectronVarManager::kR,                             3.0,  90.0, kFALSE);
  59   gammaV0Cuts->AddCut(AliDielectronVarManager::kPsiPair,                       0.0,   0.05, kFALSE);
  60   gammaV0Cuts->AddCut(AliDielectronVarManager::kM,                             0.0,   0.05, kFALSE);
  61   gammaV0Cuts->AddCut(AliDielectronVarManager::kArmPt,                         0.0,   0.05, kFALSE);
  62
  63   // selection or rejection of V0 tracks
  64   gammaV0Cuts->SetExcludeTracks(kTRUE);
  65
  66   // add the V0cuts directly to the track filter or to some cut group of it
  67
  68 */
  69 //                                                                       //
  70 ///////////////////////////////////////////////////////////////////////////
  71
  72
  73 #include "AliDielectronV0Cuts.h"
  74 #include "AliDielectronVarManager.h"
  75 #include "AliDielectronTrackCuts.h"
  76 #include "AliDielectronPID.h"
  77 #include "AliESDv0.h"
  78 #include "AliAODv0.h"
  79
  80 ClassImp(AliDielectronV0Cuts)
  81
  82
  83 AliDielectronV0Cuts::AliDielectronV0Cuts() :
  84   AliDielectronVarCuts(),
  85   fV0TrackArr(0),
  86   fExcludeTracks(kTRUE),
  87   fV0finder(kAll),
  88   fMotherPdg(0),
  89   fNegPdg(0),
  90   fPosPdg(0),
  91   fPID(-1),
  92   fPIDCutType(kBoth),
  93   fOrbit(0),
  94   fPeriod(0),
  95   fBunchCross(0)
  96 {
  97   //
  98   // Default costructor
  99   //
 100 }
 101
 102 //________________________________________________________________________
 103 AliDielectronV0Cuts::AliDielectronV0Cuts(const char* name, const char* title) :
 104   AliDielectronVarCuts(name,title),
 105   fV0TrackArr(0),
 106   fExcludeTracks(kTRUE),
 107   fV0finder(kAll),
 108   fMotherPdg(0),
 109   fNegPdg(0),
 110   fPosPdg(0),
 111   fPID(-1),
 112   fPIDCutType(kBoth),
 113   fOrbit(0),
 114   fPeriod(0),
 115   fBunchCross(0)
 116 {
 117   //
 118   // Named contructor
 119   //
 120 }
 121
 122 //________________________________________________________________________
 123 AliDielectronV0Cuts::~AliDielectronV0Cuts()
 124 {
 125   //
 126   // Destructor
 127   //
 128
 129 }
 130
 131 //________________________________________________________________________
 132 void AliDielectronV0Cuts::InitEvent(AliVTrack *trk)
 133 {
 134   //
 135   // Init the V0 candidates
 136   //
 137
 138   // take current event from the track
 139   // TODO: this should be simplyfied by AliVTrack::GetEvent() as soon as implemented
 140   const AliVEvent *ev=0;
 141   if(trk->IsA() == AliAODTrack::Class())
 142     ev=static_cast<const AliVEvent*>((static_cast<const AliAODTrack*>(trk))->GetAODEvent());
 143   else if(trk->IsA() == AliESDtrack::Class())
 144     ev=static_cast<const AliVEvent*>((static_cast<const AliESDtrack*>(trk))->GetESDEvent());
 145   else
 146     return;
 147
 148
 149   // IsNewEvent
 150   if(!ev) return;
 151   if(!IsNewEvent(ev)) return;
 152   //  printf(" Build V0 candidates according to the applied cuts \n");
 153
 154   // rest booleans
 155   fV0TrackArr.ResetAllBits();
 156
 157   // basic quality cut, /*at least one*/ both of the V0 daughters has to fullfill
 158   // always update ::Print accordingly
 159   AliDielectronVarCuts dauQAcuts1;
 160   dauQAcuts1.AddCut(AliDielectronVarManager::kNclsTPC,      70.0,  160.0);
 161   dauQAcuts1.AddCut(AliDielectronVarManager::kTPCchi2Cl,     0.0,    4.0);
 162   dauQAcuts1.AddCut(AliDielectronVarManager::kKinkIndex0,            0.0);
 163   dauQAcuts1.AddCut(AliDielectronVarManager::kEta,          -0.9,    0.9);
 164   dauQAcuts1.AddCut(AliDielectronVarManager::kPt,            0.05, 100.0);
 165   AliDielectronTrackCuts dauQAcuts2;
 166   //  dauQAcuts2.SetRequireITSRefit(kTRUE);
 167   dauQAcuts2.SetRequireTPCRefit(kTRUE);
 168   AliDielectronPID dauPIDcuts;
 169   if(fPID>=0) dauPIDcuts.SetDefaults(fPID);
 170
 171   Int_t nV0s      = 0;
 172   Int_t nV0stored = 0;
 173   AliDielectronPair candidate;
 174   candidate.SetPdgCode(fMotherPdg);
 175
 176   // ESD or AOD event
 177   if(ev->IsA() == AliESDEvent::Class()) {
 178     const AliESDEvent *esdev = static_cast<const AliESDEvent*>(ev);
 179
 180     //printf("there are %d V0s in the event \n",esdev->GetNumberOfV0s());
 181     // loop over V0s
 182     for (Int_t iv=0; iv<esdev->GetNumberOfV0s(); ++iv){
 183       AliESDv0 *v = esdev->GetV0(iv);
 184       if(!v) continue;
 185
 186       // check the v0 finder
 187       if( v->GetOnFlyStatus() && fV0finder==AliDielectronV0Cuts::kOffline  ) continue;
 188       if(!v->GetOnFlyStatus() && fV0finder==AliDielectronV0Cuts::kOnTheFly ) continue;
 189
 190       // should we make use of AliESDv0Cuts::GetPdgCode() to preselect candiadtes, e.g.:
 191       // if(fMotherPdg!=v->GetPdgCode()) continue;
 192
 193       AliESDtrack *trNeg=esdev->GetTrack(v->GetIndex(0));
 194       AliESDtrack *trPos=esdev->GetTrack(v->GetIndex(1));
 195       if(!trNeg || !trPos){
 196         printf("Error: Couldn't get V0 daughter: %p - %p\n",trNeg,trPos);
 197         continue;
 198       }
 199
 200       // protection against LS v0s
 201       if(trNeg->Charge() == trPos->Charge()) continue;
 202
 203       // PID default cuts
 204       if(fPID>=0) {
 205         Bool_t selected=kFALSE;
 206         selected=dauPIDcuts.IsSelected(trNeg);
 207         if(fPIDCutType==kBoth) selected &= dauPIDcuts.IsSelected(trPos);
 208         if(fPIDCutType==kAny)  selected |= dauPIDcuts.IsSelected(trPos);
 209         if(!selected) continue;
 210       }
 211
 212       // basic track cuts
 213       if( !dauQAcuts2.IsSelected(trNeg) ) continue;
 214       if( !dauQAcuts2.IsSelected(trPos) ) continue;
 215       if( !dauQAcuts1.IsSelected(trNeg) ) continue;
 216       if( !dauQAcuts1.IsSelected(trPos) ) continue;
 217
 218       if(fMotherPdg==22) candidate.SetGammaTracks(trNeg, 11, trPos, 11);
 219       else candidate.SetTracks(trNeg, (trNeg->Charge()<0?fNegPdg:fPosPdg), trPos, (trPos->Charge()<0?fNegPdg:fPosPdg));
 220       // eventually take the external trackparam and build the KFparticles by hand (see AliESDv0::GetKFInfo)
 221       // the folowing is not needed, because the daughters where used in the v0 vertex fit (I guess)
 222       //      AliKFVertex v0vtx = *v;
 223       //      candidate.SetProductionVertex(v0vtx);
 224
 225       if(AliDielectronVarCuts::IsSelected(&candidate)) {
 226         nV0s++;
 227         //printf(" gamma found for vtx %p dau1id %d dau2id %d \n",v,trNeg->GetID(),trPos->GetID());
 228         fV0TrackArr.SetBitNumber(trNeg->GetID());
 229         fV0TrackArr.SetBitNumber(trPos->GetID());
 230       }
 231     }
 232
 233   }
 234   else if(ev->IsA() == AliAODEvent::Class()) {
 235     const AliAODEvent *aodEv = static_cast<const AliAODEvent*>(ev);
 236     if(!aodEv->GetV0s()) return; // protection for nano AODs
 237
 238     // loop over vertices
 239     for (Int_t ivertex=0; ivertex<aodEv->GetNumberOfV0s(); ++ivertex){
 240       AliAODv0 *v=aodEv->GetV0(ivertex);
 241       if(!v) continue;
 242
 243       // check the v0 finder
 244       if( v->GetOnFlyStatus() && fV0finder==AliDielectronV0Cuts::kOffline  ) continue;
 245       if(!v->GetOnFlyStatus() && fV0finder==AliDielectronV0Cuts::kOnTheFly ) continue;
 246
 247       AliAODTrack *trNeg=dynamic_cast<AliAODTrack*>(v->GetDaughter(0));
 248       AliAODTrack *trPos=dynamic_cast<AliAODTrack*>(v->GetDaughter(1));
 249       if(!trNeg || !trPos){
 250         printf("Error: Couldn't get V0 daughter: %p - %p\n",trNeg,trPos);
 251         continue;
 252       }
 253       nV0stored++;
 254
 255       // protection against LS v0s
 256       if(trNeg->Charge() == trPos->Charge()) continue;
 257
 258       // PID default cuts
 259       if(fPID>=0) {
 260         Bool_t selected=kFALSE;
 261         selected=dauPIDcuts.IsSelected(trNeg);
 262         if(fPIDCutType==kBoth) selected &= dauPIDcuts.IsSelected(trPos);
 263         if(fPIDCutType==kAny)  selected |= dauPIDcuts.IsSelected(trPos);
 264         if(!selected) continue;
 265       }
 266
 267       // basic track cuts
 268       if( !dauQAcuts2.IsSelected(trNeg) ) continue;
 269       if( !dauQAcuts2.IsSelected(trPos) ) continue;
 270       if( !dauQAcuts1.IsSelected(trNeg) ) continue;
 271       if( !dauQAcuts1.IsSelected(trPos) ) continue;
 272
 273       AliKFVertex v0vtx = *(v->GetSecondaryVtx());
 274       if(fMotherPdg==22) candidate.SetGammaTracks(trNeg, 11, trPos, 11);
 275       else candidate.SetTracks(trNeg, (trNeg->Charge()<0?fNegPdg:fPosPdg), trPos, (trPos->Charge()<0?fNegPdg:fPosPdg));
 276       candidate.SetProductionVertex(v0vtx);
 277
 278       if(AliDielectronVarCuts::IsSelected(&candidate)) {
 279         nV0s++;
 280         //printf(" gamma found for vtx %p dau1id %d dau2id %d \n",v,trNeg->GetID(),trPos->GetID());
 281         fV0TrackArr.SetBitNumber(trNeg->GetID());
 282         fV0TrackArr.SetBitNumber(trPos->GetID());
 283       }
 284     }
 285     //    printf("there are %d V0s in the event \n",nV0stored);
 286   }
 287   else
 288     return;
 289
 290   //  printf(" Number of V0s candiates found %d/%d \n",nV0s,nV0stored);
 291
 292 }
 293 //________________________________________________________________________
 294 Bool_t AliDielectronV0Cuts::IsSelected(TObject* track)
 295 {
 296   //
 297   // Make cut decision
 298   //
 299
 300   if(!track) return kFALSE;
 301
 302   AliVTrack *vtrack = static_cast<AliVTrack*>(track);
 303   InitEvent(vtrack);
 304   //printf(" track ID %d selected result %d %d \n",vtrack->GetID(),(fV0TrackArr.TestBitNumber(vtrack->GetID())),fExcludeTracks);
 305   return ( (fV0TrackArr.TestBitNumber(vtrack->GetID()))^fExcludeTracks );
 306 }
 307
 308 //________________________________________________________________________
 309 Bool_t AliDielectronV0Cuts::IsNewEvent(const AliVEvent *ev)
 310 {
 311   //
 312   // check weather we process a new event
 313   //
 314
 315   //  printf(" current ev %d %d %d \n",fBunchCross, fOrbit, fPeriod);
 316   //  printf(" new event %p %d %d %d \n",ev, ev->GetBunchCrossNumber(), ev->GetOrbitNumber(), ev->GetPeriodNumber());
 317
 318   if( fBunchCross == ev->GetBunchCrossNumber() ) {
 319     if( fOrbit == ev->GetOrbitNumber() )         {
 320       if( fPeriod == ev->GetPeriodNumber() )     {
 321         return kFALSE;
 322       }
 323     }
 324   }
 325
 326   fBunchCross = ev->GetBunchCrossNumber();
 327   fOrbit      = ev->GetOrbitNumber();
 328   fPeriod     = ev->GetPeriodNumber();
 329   return kTRUE;
 330 }
 331
 332 //________________________________________________________________________
 333 void AliDielectronV0Cuts::Print(const Option_t* /*option*/) const
 334 {
 335   //
 336   // Print cuts and the range
 337   //
 338   printf(" V0 cuts:\n");
 339   printf(" V0 finder mode: %s \n",(fV0finder ==kOnTheFly ? "One-The-Fly":
 340                                    (fV0finder==kOffline  ? "Offline":
 341                                     "One-The-Fly+Offline") ) );
 342   AliDielectronVarCuts::Print();
 343
 344   printf(" V0 daughter cuts (applied to both):\n");
 345   AliDielectronVarCuts dauQAcuts1;
 346   dauQAcuts1.AddCut(AliDielectronVarManager::kNclsTPC,      70.0,  160.0);
 347   dauQAcuts1.AddCut(AliDielectronVarManager::kTPCchi2Cl,     0.0,    4.0);
 348   dauQAcuts1.AddCut(AliDielectronVarManager::kKinkIndex0,            0.0);
 349   dauQAcuts1.AddCut(AliDielectronVarManager::kEta,          -0.9,    0.9);
 350   dauQAcuts1.AddCut(AliDielectronVarManager::kPt,            0.05, 100.0);
 351   dauQAcuts1.Print();
 352   AliDielectronTrackCuts dauQAcuts2;
 353   //  dauQAcuts2.SetRequireITSRefit(kTRUE);
 354   dauQAcuts2.SetRequireTPCRefit(kTRUE);
 355   //  dauQAcuts2.Print(); //TODO activate as soon as implemented
 356
 357   if(fPID>=0) {
 358     printf(" V0 daughter PID cuts (applied to %s):\n",(fPIDCutType==kBoth?"both":"any"));
 359     AliDielectronPID dauPIDcuts;
 360     dauPIDcuts.SetDefaults(fPID);
 361     dauPIDcuts.Print(); //TODO activate as soon as implemented
 362   }
 363
 364 }