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 (fTrackEfficiency < 1) {
 133           Double_t r = gRandom->Rndm();
 134           if (fTrackEfficiency < r)
 135             continue;
 136         }
 137
 138         if (!fEsdTrackCuts->AcceptTrack(etrack))
 139           continue;
 140
 141         AliESDtrack *ntrack = AliESDtrackCuts::GetTPCOnlyTrack(fEsdEv,etrack->GetID());
 142         if (!ntrack)
 143           continue;
 144         if (ntrack->Pt()<=0) {
 145           delete ntrack;
 146           continue;
 147         }
 148         Double_t bfield[3] = {0,0,0};
 149         ntrack->GetBxByBz(bfield);
 150         AliExternalTrackParam exParam;
 151         Bool_t relate = ntrack->RelateToVertexBxByBz(vtxSPD,bfield,kVeryBig,&exParam);
 152         if (!relate) {
 153           delete ntrack;
 154           continue;
 155         }
 156         // set the constraint parameters to the track
 157         ntrack->Set(exParam.GetX(),exParam.GetAlpha(),exParam.GetParameter(),exParam.GetCovariance());
 158         if (ntrack->Pt()<=0) {
 159           delete ntrack;
 160           continue;
 161         }
 162         if (fDoPropagation)
 163           AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(ntrack,fDist);
 164         new ((*fTracks)[ntrnew++]) AliESDtrack(*ntrack);
 165         delete ntrack;
 166       }
 167     } else { /* no spd vtx constraint */
 168       Int_t ntr = fEsdEv->GetNumberOfTracks();
 169       for (Int_t i=0, ntrnew=0; i<ntr; ++i) {
 170         AliESDtrack *etrack = fEsdEv->GetTrack(i);
 171         if (!etrack)
 172           continue;
 173         if (fTrackEfficiency < 1) {
 174           Double_t r = gRandom->Rndm();
 175           if (fTrackEfficiency < r)
 176             continue;
 177         }
 178         if ((fEsdTrackCuts!=0) && !fEsdTrackCuts->AcceptTrack(etrack))
 179           continue;
 180         AliESDtrack *ntrack = new ((*fTracks)[ntrnew++]) AliESDtrack(*etrack);
 181         if (fDoPropagation)
 182           AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(ntrack,fDist);
 183       }
 184     }
 185
 186   } else { // use hybrid track cuts
 187
 188     am->LoadBranch("Tracks");
 189     Int_t ntr = fEsdEv->GetNumberOfTracks();
 190     for (Int_t i=0, ntrnew=0; i<ntr; ++i) {
 191       AliESDtrack *etrack = fEsdEv->GetTrack(i);
 192       if (!etrack)
 193         continue;
 194
 195       if (fTrackEfficiency < 1) {
 196         Double_t r = gRandom->Rndm();
 197         if (fTrackEfficiency < r)
 198           continue;
 199       }
 200
 201       if (fEsdTrackCuts->AcceptTrack(etrack)) {
 202         AliESDtrack *newTrack = new ((*fTracks)[ntrnew]) AliESDtrack(*etrack);
 203         if (fDoPropagation)
 204           AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(newTrack,fDist);
 205         newTrack->SetBit(BIT(22),0);
 206         newTrack->SetBit(BIT(23),0);
 207         ++ntrnew;
 208       } else if (fHybridTrackCuts->AcceptTrack(etrack)) {
 209         if (!etrack->GetConstrainedParam())
 210           continue;
 211         UInt_t status = etrack->GetStatus();
 212         if (!fIncludeNoITS && ((status&AliESDtrack::kITSrefit)==0))
 213           continue;
 214         AliESDtrack *newTrack = new ((*fTracks)[ntrnew]) AliESDtrack(*etrack);
 215         if (fDoPropagation)
 216           AliEMCALRecoUtils::ExtrapolateTrackToEMCalSurface(newTrack,fDist);
 217         const AliExternalTrackParam* constrainParam = etrack->GetConstrainedParam();
 218         newTrack->Set(constrainParam->GetX(),
 219                       constrainParam->GetAlpha(),
 220                       constrainParam->GetParameter(),
 221                       constrainParam->GetCovariance());
 222         if ((status&AliESDtrack::kITSrefit)==0) {
 223           newTrack->SetBit(BIT(22),0); //type 2
 224           newTrack->SetBit(BIT(23),1);
 225         } else {
 226           newTrack->SetBit(BIT(22),1); //type 1
 227           newTrack->SetBit(BIT(23),0);
 228         }
 229         ++ntrnew;
 230       }
 231     }
 232   }
 233 }