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