PWG4/UserTasks/EmcalTasks/AliEmcalEsdTpcTrackTask.cxx

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