PWG/EMCAL/AliEmcalEsdTrackFilterTask.cxx

   1 //
   2 // Task to filter Esd tracks and propagate to Emcal surface.
   3 //
   4 // Author: C.Loizides
   5
   6 #include "AliEmcalEsdTrackFilterTask.h"
   7 #include <TClonesArray.h>
   8 #include <TRandom3.h>
   9 #include <TGeoGlobalMagField.h>
  10 #include <AliAnalysisManager.h>
  11 #include <AliEMCALRecoUtils.h>
  12 #include <AliESDEvent.h>
  13 #include <AliESDtrackCuts.h>
  14 #include <AliMagF.h>
  15 #include <AliTrackerBase.h>
  16
  17
  18 ClassImp(AliEmcalEsdTrackFilterTask)
  19
  20 //________________________________________________________________________
  21 AliEmcalEsdTrackFilterTask::AliEmcalEsdTrackFilterTask() :
  22   AliAnalysisTaskSE("AliEmcalEsdTrackFilterTask"),
  23   fEsdTrackCuts(0),
  24   fDoSpdVtxCon(0),
  25   fHybridTrackCuts(0),
  26   fTracksName(),
  27   fIncludeNoITS(kTRUE),
  28   fDoPropagation(kFALSE),
  29   fDist(440),
  30   fTrackEfficiency(0),
  31   fIsMC(kFALSE),
  32   fEsdEv(0),
  33   fTracks(0)
  34 {
  35   // Constructor.
  36 }
  37
  38 //________________________________________________________________________
  39 AliEmcalEsdTrackFilterTask::AliEmcalEsdTrackFilterTask(const char *name) :
  40   AliAnalysisTaskSE(name),
  41   fEsdTrackCuts(0),
  42   fDoSpdVtxCon(0),
  43   fHybridTrackCuts(0),
  44   fTracksName("EsdTracksOut"),
  45   fIncludeNoITS(kTRUE),
  46   fDoPropagation(kFALSE),
  47   fDist(440),
  48   fTrackEfficiency(0),
  49   fIsMC(kFALSE),
  50   fEsdEv(0),
  51   fTracks(0)
  52 {
  53   // Constructor.
  54
  55   if (!name)
  56     return;
  57
  58   SetName(name);
  59
  60   fBranchNames = "ESD:AliESDHeader.,AliESDRun.,SPDVertex.,Tracks";
  61 }
  62
  63 //________________________________________________________________________
  64 AliEmcalEsdTrackFilterTask::~AliEmcalEsdTrackFilterTask()
  65 {
  66   //Destructor
  67
  68   delete fEsdTrackCuts;
  69 }
  70
  71 //________________________________________________________________________
  72 void AliEmcalEsdTrackFilterTask::UserCreateOutputObjects()
  73 {
  74   // Create histograms.
  75
  76   fTracks = new TClonesArray("AliESDtrack");
  77   fTracks->SetName(fTracksName);
  78
  79   if (fDoSpdVtxCon) {
  80     if (!fEsdTrackCuts) {
  81       AliInfo("No track cuts given, creating default (standard only TPC) cuts");
  82       fEsdTrackCuts = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
  83       fEsdTrackCuts->SetPtRange(0.15,1e3);
  84     }
  85   } else {
  86     AliWarning("No track cuts given, but maybe this is indeed intended?");
  87   }
  88 }
  89
  90 //________________________________________________________________________
  91 void AliEmcalEsdTrackFilterTask::UserExec(Option_t *)
  92 {
  93   // Main loop, called for each event.
  94
  95   fEsdEv = dynamic_cast<AliESDEvent*>(InputEvent());
  96   if (!fEsdEv) {
  97     AliError("Task works only on ESD events, returning");
  98     return;
  99   }
 100
 101   AliAnalysisManager *am = AliAnalysisManager::GetAnalysisManager();
 102   if (!am) {
 103     AliError("Manager zero, returning");
 104     return;
 105   }
 106
 107   // add tracks to event if not yet there
 108   fTracks->Delete();
 109   if (!(InputEvent()->FindListObject(fTracksName)))
 110     InputEvent()->AddObject(fTracks);
 111
 112   if (!fHybridTrackCuts) { // constrain TPC tracks to SPD vertex if fDoSpdVtxCon==kTRUE
 113     am->LoadBranch("AliESDRun.");
 114     am->LoadBranch("AliESDHeader.");
 115     am->LoadBranch("Tracks");
 116
 117     if (fDoSpdVtxCon) {
 118       if (!TGeoGlobalMagField::Instance()->GetField()) { // construct field map
 119         fEsdEv->InitMagneticField();
 120       }
 121       am->LoadBranch("SPDVertex.");
 122       const AliESDVertex *vtxSPD = fEsdEv->GetPrimaryVertexSPD();
 123       if (!vtxSPD) {
 124         AliError("No SPD vertex, returning");
 125         return;
 126       }
 127       Int_t ntr = fEsdEv->GetNumberOfTracks();
 128       for (Int_t i=0, ntrnew=0; i<ntr; ++i) {
 129         AliESDtrack *etrack = fEsdEv->GetTrack(i);
 130         if (!etrack)
 131           continue;
 132
 133         if (!fEsdTrackCuts->AcceptTrack(etrack))
 134           continue;
 135
 136         AliESDtrack *ntrack = AliESDtrackCuts::GetTPCOnlyTrack(fEsdEv,etrack->GetID());
 137         if (!ntrack)
 138           continue;
 139         if (ntrack->Pt()<=0) {
 140           delete ntrack;
 141           continue;
 142         }
 143         Double_t bfield[3] = {0,0,0};
 144         ntrack->GetBxByBz(bfield);
 145         AliExternalTrackParam exParam;
 146         Bool_t relate = ntrack->RelateToVertexBxByBz(vtxSPD,bfield,kVeryBig,&exParam);
 147         if (!relate) {
 148           delete ntrack;
 149           continue;
 150         }
 151         // set the constraint parameters to the track
 152         ntrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 153         if (ntrack->Pt()<=0) {
 154           delete ntrack;
 155           continue;
 156         }
 157
 158         if (fTrackEfficiency) {
 159           Double_t r = gRandom->Rndm();
 160           if (fTrackEfficiency->Eval(ntrack->Pt()) < r)
 161             continue;
 162         }
 163
 164         if (fDoPropagation)
 165           AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(ntrack,fDist);
 166         new ((*fTracks)[ntrnew++]) AliESDtrack(*ntrack);
 167         delete ntrack;
 168       }
 169     } else { /* no spd vtx constraint */
 170       Int_t ntr = fEsdEv->GetNumberOfTracks();
 171       for (Int_t i=0, ntrnew=0; i<ntr; ++i) {
 172         AliESDtrack *etrack = fEsdEv->GetTrack(i);
 173         if (!etrack)
 174           continue;
 175
 176         if ((fEsdTrackCuts!=0) && !fEsdTrackCuts->AcceptTrack(etrack))
 177           continue;
 178
 179         if (fTrackEfficiency) {
 180           Double_t r = gRandom->Rndm();
 181           if (fTrackEfficiency->Eval(etrack->Pt()) < r)
 182             continue;
 183         }
 184
 185         AliESDtrack *ntrack = new ((*fTracks)[ntrnew++]) AliESDtrack(*etrack);
 186         if (fDoPropagation)
 187           AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(ntrack,fDist);
 188       }
 189     }
 190
 191   } else { // use hybrid track cuts
 192
 193     am->LoadBranch("Tracks");
 194     Int_t ntr = fEsdEv->GetNumberOfTracks();
 195     for (Int_t i=0, ntrnew=0; i<ntr; ++i) {
 196       AliESDtrack *etrack = fEsdEv->GetTrack(i);
 197       if (!etrack)
 198         continue;
 199
 200       if (fEsdTrackCuts->AcceptTrack(etrack)) {
 201         if (fTrackEfficiency) {
 202           Double_t r = gRandom->Rndm();
 203           if (fTrackEfficiency->Eval(etrack->Pt()) < r)
 204             continue;
 205         }
 206         AliESDtrack *newTrack = new ((*fTracks)[ntrnew]) AliESDtrack(*etrack);
 207         if (fDoPropagation)
 208           AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(newTrack,fDist);
 209         newTrack->SetBit(BIT(22),0);
 210         newTrack->SetBit(BIT(23),0);
 211         if (!fMCEvent) newTrack->SetLabel(0);
 212         ++ntrnew;
 213       } else if (fHybridTrackCuts->AcceptTrack(etrack)) {
 214         if (!etrack->GetConstrainedParam())
 215           continue;
 216         UInt_t status = etrack->GetStatus();
 217         if (!fIncludeNoITS && ((status&AliESDtrack::kITSrefit)==0))
 218           continue;
 219
 220         if (fTrackEfficiency) {
 221           Double_t r = gRandom->Rndm();
 222           if (fTrackEfficiency->Eval(etrack->Pt()) < r)
 223             continue;
 224         }
 225         AliESDtrack *newTrack = new ((*fTracks)[ntrnew]) AliESDtrack(*etrack);
 226         if (fDoPropagation)
 227           AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(newTrack,fDist);
 228         const AliExternalTrackParam* constrainParam = etrack->GetConstrainedParam();
 229         newTrack->Set(constrainParam->GetX(),
 230                       constrainParam->GetAlpha(),
 231                       constrainParam->GetParameter(),
 232                       constrainParam->GetCovariance());
 233         if ((status&AliESDtrack::kITSrefit)==0) {
 234           newTrack->SetBit(BIT(22),0); //type 2
 235           newTrack->SetBit(BIT(23),1);
 236         } else {
 237           newTrack->SetBit(BIT(22),1); //type 1
 238           newTrack->SetBit(BIT(23),0);
 239         }
 240         if (!fMCEvent) newTrack->SetLabel(0);
 241         ++ntrnew;
 242       }
 243     }
 244   }
 245 }