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