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, fMCEvent);
 156   if (!fEventCuts.IsSelected(&event)) return;
 157
 158   // check the event
 159   Int_t type = EventEval(esd);
 160
 161   // if processable, then process it
 162   if      (type == 0) ProcessESD(esd, esd->GetPrimaryVertexTracks(), stack);
 163   else if (type == 1) ProcessESD(esd, esd->GetPrimaryVertexSPD()   , stack);
 164   else return;
 165
 166   // update histogram container
 167   PostData(1, fOut);
 168 }
 169
 170 //__________________________________________________________________________________________________
 171 void AliRsnAnalysisMonitorTask::Terminate(Option_t *)
 172 {
 173 //
 174 // Terminate
 175 //
 176 }
 177
 178 //__________________________________________________________________________________________________
 179 Int_t AliRsnAnalysisMonitorTask::EventEval(AliESDEvent *esd)
 180 {
 181 //
 182 // Checks if the event is good for analysis.
 183 // Returns:
 184 // ---> 0 if a good primary vertex with tracks was found,
 185 // ---> 1 if a good SPD primary vertex was found
 186 // ---> 2 otherwise (event to be rejected)
 187 // In any case, adds an entry to the TTree, to keep trace of all events.
 188 //
 189
 190   // get the best primary vertex:
 191   // first try that with tracks, then the SPD one
 192   const AliESDVertex *vTrk  = esd->GetPrimaryVertexTracks();
 193   const AliESDVertex *vSPD  = esd->GetPrimaryVertexSPD();
 194   if(vTrk->GetNContributors() > 0)
 195   {
 196     return 0;
 197   }
 198   else if (vSPD->GetNContributors() > 0)
 199   {
 200     return 1;
 201   }
 202   else
 203   {
 204     return 2;
 205   }
 206 }
 207
 208 //__________________________________________________________________________________________________
 209 void AliRsnAnalysisMonitorTask::ProcessESD
 210 (AliESDEvent *esd, const AliESDVertex *v, AliStack *stack)
 211 {
 212 //
 213 // Process the ESD container, to read all tracks and copy their useful values.
 214 // All info are stored into an AliRsnMonitorTrack object and saved into the
 215 // TClonesArray which is one of the branches of the output TTree.
 216 //
 217
 218   // create interfacr objects
 219   AliRsnEvent    event;
 220   AliRsnDaughter daughter;
 221   event.SetRef(esd, fMCEvent);
 222
 223   // ITS stuff #1
 224   // create the response function and initialize it to MC or not
 225   // depending if the AliStack object is there or not
 226   Bool_t isMC = (stack != 0x0);
 227   AliITSPIDResponse itsrsp(isMC);
 228
 229   // TOF stuff #1: init OCDB
 230   Int_t run = esd->GetRunNumber();
 231   AliCDBManager *cdb = AliCDBManager::Instance();
 232   cdb->SetDefaultStorage("raw://");
 233   cdb->SetRun(run);
 234   // TOF stuff #2: init calibration
 235   fTOFcalib->SetCorrectTExp(fTOFcorrectTExp);
 236   fTOFcalib->Init();
 237   // TOF stuff #3: calibrate
 238   if (fTOFcalibrateESD) fTOFcalib->CalibrateESD(esd);
 239   if (fTOFtuneMC) fTOFmaker->TuneForMC(esd);
 240   if (fTOFuseT0)
 241   {
 242     fTOFmaker->ComputeT0TOF(esd);
 243     fTOFmaker->ApplyT0TOF(esd);
 244     fESDpid->MakePID(esd, kFALSE, 0.);
 245   }
 246
 247   // loop on all tracks
 248   Int_t               i, k, nITS, ntracks = esd->GetNumberOfTracks();;
 249   Bool_t              isTPC, isITSSA, isTOF;
 250   Float_t             b[2], bCov[3];
 251   Double_t            time[10];
 252
 253   for (i = 0; i < ntracks; i++)
 254   {
 255     AliESDtrack *track = esd->GetTrack(i);
 256     event.SetDaughter(daughter, i, AliRsnDaughter::kTrack);
 257
 258     // reset the output object
 259     // 'usable' flag will need to be set to 'ok'
 260     fTrack->Reset();
 261
 262     // check cuts
 263     fTrack->CutsPassed() = fTrackCuts.IsSelected(&daughter);
 264
 265     // skip NULL pointers, kink daughters and tracks which
 266     // cannot be propagated to primary vertex
 267     if (!track) continue;
 268     if ((Int_t)track->GetKinkIndex(0) > 0) continue;
 269     if (!track->RelateToVertex(v, esd->GetMagneticField(), kVeryBig)) continue;
 270
 271     // get MC info if possible
 272     if (stack) fTrack->AdoptMC(TMath::Abs(track->GetLabel()), stack);
 273
 274     // copy general info
 275     fTrack->Status() = (UInt_t)track->GetStatus();
 276     fTrack->Length() = (Double_t)track->GetIntegratedLength();
 277     fTrack->Charge() = (Int_t)track->Charge();
 278     fTrack->PrecX()  = (Double_t)track->Px();
 279     fTrack->PrecY()  = (Double_t)track->Py();
 280     fTrack->PrecZ()  = (Double_t)track->Pz();
 281
 282     // evaluate some flags from the status to decide what to do next in some points
 283     isTPC   = ((fTrack->Status() & AliESDtrack::kTPCin)  != 0);
 284     isITSSA = ((fTrack->Status() & AliESDtrack::kTPCin)  == 0 && (fTrack->Status() & AliESDtrack::kITSrefit) != 0 && (fTrack->Status() & AliESDtrack::kITSpureSA) == 0 && (fTrack->Status() & AliESDtrack::kITSpid) != 0);
 285     isTOF   = ((fTrack->Status() & AliESDtrack::kTOFout) != 0 && (fTrack->Status() & AliESDtrack::kTIME) != 0);
 286
 287     // accept only tracks which are TPC+ITS or ITS standalone
 288     if (isITSSA)
 289     {
 290       fTrack->ITSsa() = kTRUE;
 291       fTrack->TOFok() = kFALSE;
 292     }
 293     else if (isTPC)
 294     {
 295       fTrack->ITSsa() = kFALSE;
 296       fTrack->TOFok() = isTOF;
 297     }
 298     else
 299       continue;
 300
 301     // get DCA to primary vertex
 302     track->GetImpactParameters(b, bCov);
 303     fTrack->DCAr() = (Double_t)b[0];
 304     fTrack->DCAz() = (Double_t)b[1];
 305
 306     // get ITS info
 307     for (k = 0; k < 6; k++)
 308     {
 309       fTrack->ITSmap(k) = track->HasPointOnITSLayer(k);
 310     }
 311     if (isITSSA)
 312     {
 313       fTrack->ITSchi2() = track->GetITSchi2();
 314       fTrack->ITSsignal() = track->GetITSsignal();
 315       nITS = fTrack->SSDcount() + fTrack->SDDcount();
 316       for (k = 0; k < AliPID::kSPECIES; k++)
 317       {
 318         fTrack->ITSnsigma(k) = itsrsp.GetNumberOfSigmas(fTrack->Prec(), fTrack->ITSsignal(), (AliPID::EParticleType)k, nITS, kTRUE);
 319       }
 320     }
 321
 322     // get TPC info
 323     if (isTPC)
 324     {
 325       fTrack->TPCcount()  = (Int_t)track->GetTPCclusters(0);
 326       fTrack->TPCchi2()   = (Double_t)track->GetTPCchi2();
 327       fTrack->TPCsignal() = (Double_t)track->GetTPCsignal();
 328       fTrack->PtpcX()     = fTrack->PtpcY() = fTrack->PtpcZ() = 1E10;
 329       if (track->GetInnerParam())
 330       {
 331         fTrack->PtpcX() = track->GetInnerParam()->Px();
 332         fTrack->PtpcY() = track->GetInnerParam()->Py();
 333         fTrack->PtpcZ() = track->GetInnerParam()->Pz();
 334         for (k = 0; k < AliPID::kSPECIES; k++)
 335         {
 336           fTrack->TPCnsigma(k) = fESDpid->NumberOfSigmasTPC(track, (AliPID::EParticleType)k);
 337         }
 338       }
 339     }
 340
 341     // get TOF info
 342     if (isTOF)
 343     {
 344       track->GetIntegratedTimes(time);
 345       fTrack->TOFsignal() = (Double_t)track->GetTOFsignal();
 346       for (k = 0; k < AliPID::kSPECIES; k++)
 347       {
 348         fTrack->TOFref(k)   = time[k];
 349         fTrack->TOFsigma(k) = (Double_t)fTOFmaker->GetExpectedSigma(fTrack->Prec(), time[k], AliPID::ParticleMass(k));
 350       }
 351     }
 352
 353     // add entry to TTree
 354     fOut->Fill();
 355   }
 356 }