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