PWG2/RESONANCES/AliRsnAnalysisMonitorTask.cxx

   1 //
   2 // Implementation file for implementation of data analysis aft 900 GeV
   3 //
   4 // Author: A. Pulvirenti
   5 //
   6
   7 #include "Riostream.h"
   8 #include <iomanip>
   9
  10 #include "TH1.h"
  11 #include "TTree.h"
  12 #include "TParticle.h"
  13 #include "TRandom.h"
  14 #include "TLorentzVector.h"
  15
  16 #include "AliLog.h"
  17 #include "AliESDpid.h"
  18 #include "AliESDEvent.h"
  19 #include "AliESDVertex.h"
  20 #include "AliESDtrack.h"
  21 #include "AliStack.h"
  22 #include "AliMCEvent.h"
  23 #include "AliTOFT0maker.h"
  24 #include "AliTOFcalib.h"
  25 #include "AliCDBManager.h"
  26 #include "AliITSPIDResponse.h"
  27 #include "AliRsnMonitorTrack.h"
  28 #include "AliRsnDaughter.h"
  29 #include "AliRsnEvent.h"
  30
  31 #include "AliRsnAnalysisMonitorTask.h"
  32
  33 ClassImp(AliRsnAnalysisMonitorTask)
  34
  35 //__________________________________________________________________________________________________
  36 AliRsnAnalysisMonitorTask::AliRsnAnalysisMonitorTask(const char *name) :
  37    AliAnalysisTaskSE(name),
  38    fOut(0x0),
  39    fTrack(0x0),
  40    fESDpid(0x0),
  41    fTOFmaker(0x0),
  42    fTOFcalib(0x0),
  43    fTOFcalibrateESD(kFALSE),
  44    fTOFcorrectTExp(kFALSE),
  45    fTOFuseT0(kFALSE),
  46    fTOFtuneMC(kFALSE),
  47    fTOFresolution(0.0),
  48    fEventCuts("eventCuts", AliRsnCut::kEvent),
  49    fTrackCuts("trackCuts", AliRsnCut::kDaughter)
  50
  51 {
  52 //
  53 // Constructor
  54 //
  55
  56    DefineOutput(1, TTree::Class());
  57 }
  58
  59 //__________________________________________________________________________________________________
  60 AliRsnAnalysisMonitorTask::AliRsnAnalysisMonitorTask(const AliRsnAnalysisMonitorTask& copy) :
  61    AliAnalysisTaskSE(copy),
  62    fOut(0x0),
  63    fTrack(0x0),
  64    fESDpid(0x0),
  65    fTOFmaker(0x0),
  66    fTOFcalib(0x0),
  67    fTOFcalibrateESD(kFALSE),
  68    fTOFcorrectTExp(kFALSE),
  69    fTOFuseT0(kFALSE),
  70    fTOFtuneMC(kFALSE),
  71    fTOFresolution(0.0),
  72    fEventCuts(copy.fEventCuts),
  73    fTrackCuts(copy.fTrackCuts)
  74 {
  75 //
  76 // Copy constructor
  77 //
  78 }
  79
  80 //__________________________________________________________________________________________________
  81 AliRsnAnalysisMonitorTask& AliRsnAnalysisMonitorTask::operator=(const AliRsnAnalysisMonitorTask& copy)
  82 {
  83 //
  84 // Assignment operator
  85 //
  86
  87    fTOFcalibrateESD = copy.fTOFcalibrateESD;
  88    fTOFcorrectTExp = copy.fTOFcorrectTExp;
  89    fTOFuseT0 = copy.fTOFuseT0;
  90    fTOFtuneMC = copy.fTOFtuneMC;
  91    fTOFresolution = copy.fTOFresolution;
  92
  93    return (*this);
  94 }
  95
  96 //__________________________________________________________________________________________________
  97 AliRsnAnalysisMonitorTask::~AliRsnAnalysisMonitorTask()
  98 {
  99 //
 100 // Destructor
 101 //
 102
 103    if (fOut)    delete fOut;
 104    if (fESDpid) delete fESDpid;
 105    if (fTrack)  delete fTrack;
 106 }
 107
 108 //__________________________________________________________________________________________________
 109 void AliRsnAnalysisMonitorTask::UserCreateOutputObjects()
 110 {
 111 //
 112 // Create the output data container
 113 //
 114
 115    // setup TPC response
 116    fESDpid = new AliESDpid;
 117    fESDpid->GetTPCResponse().SetBetheBlochParameters(fTPCpar[0], fTPCpar[1], fTPCpar[2], fTPCpar[3], fTPCpar[4]);
 118
 119    // setup TOF maker & calibration
 120    fTOFcalib = new AliTOFcalib;
 121    fTOFmaker = new AliTOFT0maker(fESDpid, fTOFcalib);
 122    fTOFmaker->SetTimeResolution(fTOFresolution);
 123
 124    // create output branch object
 125    fTrack = new AliRsnMonitorTrack;
 126
 127    // create output tree
 128    OpenFile(1);
 129    fOut = new TTree("rsnMonitor", "Informations on single tracks for cut checking");
 130    fOut->Branch("tracks", "AliRsnMonitorTrack", &fTrack);
 131 }
 132
 133 //__________________________________________________________________________________________________
 134 void AliRsnAnalysisMonitorTask::UserExec(Option_t *)
 135 {
 136 //
 137 // Main execution function.
 138 // Fills the fHEvents data member with the following legenda:
 139 // 0 -- event OK, prim vertex with tracks
 140 // 1 -- event OK, prim vertex with SPD
 141 // 2 -- event OK but vz large
 142 // 3 -- event bad
 143 //
 144
 145    // retrieve ESD event and related stack (if available)
 146    AliESDEvent *esd   = dynamic_cast<AliESDEvent*>(fInputEvent);
 147    AliStack    *stack = 0x0;
 148
 149    // skip NULL events
 150    if (!esd) return;
 151    if (fMCEvent) stack = fMCEvent->Stack();
 152
 153    // create interface objects to AliRsnEvent to check event cuts
 154    AliRsnEvent event;
 155    event.SetRef(esd);
 156    event.SetRefMC(fMCEvent);
 157    if (!fEventCuts.IsSelected(&event)) return;
 158
 159    // check the event
 160    Int_t type = EventEval(esd);
 161
 162    // if processable, then process it
 163    if (type == 0) ProcessESD(esd, esd->GetPrimaryVertexTracks(), stack);
 164    else if (type == 1) ProcessESD(esd, esd->GetPrimaryVertexSPD()   , stack);
 165    else return;
 166
 167    // update histogram container
 168    PostData(1, fOut);
 169 }
 170
 171 //__________________________________________________________________________________________________
 172 void AliRsnAnalysisMonitorTask::Terminate(Option_t *)
 173 {
 174 //
 175 // Terminate
 176 //
 177 }
 178
 179 //__________________________________________________________________________________________________
 180 Int_t AliRsnAnalysisMonitorTask::EventEval(AliESDEvent *esd)
 181 {
 182 //
 183 // Checks if the event is good for analysis.
 184 // Returns:
 185 // ---> 0 if a good primary vertex with tracks was found,
 186 // ---> 1 if a good SPD primary vertex was found
 187 // ---> 2 otherwise (event to be rejected)
 188 // In any case, adds an entry to the TTree, to keep trace of all events.
 189 //
 190
 191    // get the best primary vertex:
 192    // first try that with tracks, then the SPD one
 193    const AliESDVertex *vTrk  = esd->GetPrimaryVertexTracks();
 194    const AliESDVertex *vSPD  = esd->GetPrimaryVertexSPD();
 195    if (vTrk->GetNContributors() > 0) {
 196       return 0;
 197    } else if (vSPD->GetNContributors() > 0) {
 198       return 1;
 199    } else {
 200       return 2;
 201    }
 202 }
 203
 204 //__________________________________________________________________________________________________
 205 void AliRsnAnalysisMonitorTask::ProcessESD
 206 (AliESDEvent *esd, const AliESDVertex *v, AliStack *stack)
 207 {
 208 //
 209 // Process the ESD container, to read all tracks and copy their useful values.
 210 // All info are stored into an AliRsnMonitorTrack object and saved into the
 211 // TClonesArray which is one of the branches of the output TTree.
 212 //
 213
 214    // create interfacr objects
 215    AliRsnEvent    event;
 216    AliRsnDaughter daughter;
 217    event.SetRef(esd);
 218    event.SetRefMC(fMCEvent);
 219
 220    // ITS stuff #1
 221    // create the response function and initialize it to MC or not
 222    // depending if the AliStack object is there or not
 223    Bool_t isMC = (stack != 0x0);
 224    AliITSPIDResponse itsrsp(isMC);
 225
 226    // TOF stuff #1: init OCDB
 227    Int_t run = esd->GetRunNumber();
 228    AliCDBManager *cdb = AliCDBManager::Instance();
 229    cdb->SetDefaultStorage("raw://");
 230    cdb->SetRun(run);
 231    // TOF stuff #2: init calibration
 232    fTOFcalib->SetCorrectTExp(fTOFcorrectTExp);
 233    fTOFcalib->Init();
 234    // TOF stuff #3: calibrate
 235    if (fTOFcalibrateESD) fTOFcalib->CalibrateESD(esd);
 236    if (fTOFtuneMC) fTOFmaker->TuneForMC(esd);
 237    if (fTOFuseT0) {
 238       fTOFmaker->ComputeT0TOF(esd);
 239       fTOFmaker->ApplyT0TOF(esd);
 240       fESDpid->MakePID(esd, kFALSE, 0.);
 241    }
 242
 243    // loop on all tracks
 244    Int_t               i, k, nITS, ntracks = esd->GetNumberOfTracks();;
 245    Bool_t              isTPC, isITSSA, isTOF;
 246    Float_t             b[2], bCov[3];
 247    Double_t            time[10];
 248
 249    for (i = 0; i < ntracks; i++) {
 250       AliESDtrack *track = esd->GetTrack(i);
 251       event.SetDaughter(daughter, i, AliRsnDaughter::kTrack);
 252
 253       // reset the output object
 254       // 'usable' flag will need to be set to 'ok'
 255       fTrack->Reset();
 256
 257       // check cuts
 258       fTrack->CutsPassed() = fTrackCuts.IsSelected(&daughter);
 259
 260       // skip NULL pointers, kink daughters and tracks which
 261       // cannot be propagated to primary vertex
 262       if (!track) continue;
 263       if ((Int_t)track->GetKinkIndex(0) > 0) continue;
 264       if (!track->RelateToVertex(v, esd->GetMagneticField(), kVeryBig)) continue;
 265
 266       // get MC info if possible
 267       if (stack) fTrack->AdoptMC(TMath::Abs(track->GetLabel()), stack);
 268
 269       // copy general info
 270       fTrack->Status() = (UInt_t)track->GetStatus();
 271       fTrack->Length() = (Double_t)track->GetIntegratedLength();
 272       fTrack->Charge() = (Int_t)track->Charge();
 273       fTrack->PrecX()  = (Double_t)track->Px();
 274       fTrack->PrecY()  = (Double_t)track->Py();
 275       fTrack->PrecZ()  = (Double_t)track->Pz();
 276
 277       // evaluate some flags from the status to decide what to do next in some points
 278       isTPC   = ((fTrack->Status() & AliESDtrack::kTPCin)  != 0);
 279       isITSSA = ((fTrack->Status() & AliESDtrack::kTPCin)  == 0 && (fTrack->Status() & AliESDtrack::kITSrefit) != 0 && (fTrack->Status() & AliESDtrack::kITSpureSA) == 0 && (fTrack->Status() & AliESDtrack::kITSpid) != 0);
 280       isTOF   = ((fTrack->Status() & AliESDtrack::kTOFout) != 0 && (fTrack->Status() & AliESDtrack::kTIME) != 0);
 281
 282       // accept only tracks which are TPC+ITS or ITS standalone
 283       if (isITSSA) {
 284          fTrack->ITSsa() = kTRUE;
 285          fTrack->TOFok() = kFALSE;
 286       } else if (isTPC) {
 287          fTrack->ITSsa() = kFALSE;
 288          fTrack->TOFok() = isTOF;
 289       } else
 290          continue;
 291
 292       // get DCA to primary vertex
 293       track->GetImpactParameters(b, bCov);
 294       fTrack->DCAr() = (Double_t)b[0];
 295       fTrack->DCAz() = (Double_t)b[1];
 296
 297       // get ITS info
 298       for (k = 0; k < 6; k++) {
 299          fTrack->ITSmap(k) = track->HasPointOnITSLayer(k);
 300       }
 301       if (isITSSA) {
 302          fTrack->ITSchi2() = track->GetITSchi2();
 303          fTrack->ITSsignal() = track->GetITSsignal();
 304          nITS = fTrack->SSDcount() + fTrack->SDDcount();
 305          for (k = 0; k < AliPID::kSPECIES; k++) {
 306             fTrack->ITSnsigma(k) = itsrsp.GetNumberOfSigmas(fTrack->Prec(), fTrack->ITSsignal(), (AliPID::EParticleType)k, nITS, kTRUE);
 307          }
 308       }
 309
 310       // get TPC info
 311       if (isTPC) {
 312          fTrack->TPCcount()  = (Int_t)track->GetTPCclusters(0);
 313          fTrack->TPCchi2()   = (Double_t)track->GetTPCchi2();
 314          fTrack->TPCsignal() = (Double_t)track->GetTPCsignal();
 315          fTrack->PtpcX()     = fTrack->PtpcY() = fTrack->PtpcZ() = 1E10;
 316          if (track->GetInnerParam()) {
 317             fTrack->PtpcX() = track->GetInnerParam()->Px();
 318             fTrack->PtpcY() = track->GetInnerParam()->Py();
 319             fTrack->PtpcZ() = track->GetInnerParam()->Pz();
 320             for (k = 0; k < AliPID::kSPECIES; k++) {
 321                fTrack->TPCnsigma(k) = fESDpid->NumberOfSigmasTPC(track, (AliPID::EParticleType)k);
 322             }
 323          }
 324       }
 325
 326       // get TOF info
 327       if (isTOF) {
 328          track->GetIntegratedTimes(time);
 329          fTrack->TOFsignal() = (Double_t)track->GetTOFsignal();
 330          for (k = 0; k < AliPID::kSPECIES; k++) {
 331             fTrack->TOFref(k)   = time[k];
 332             fTrack->TOFsigma(k) = (Double_t)fTOFmaker->GetExpectedSigma(fTrack->Prec(), time[k], AliPID::ParticleMass(k));
 333          }
 334       }
 335
 336       // add entry to TTree
 337       fOut->Fill();
 338    }
 339 }