PWG/EMCAL/AliEmcalEsdTrackFilterTask.cxx

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