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