PWG/EMCAL/AliEmcalEsdTpcTrackTask.cxx

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