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
  29 #include "AliRsnAnalysisMonitorTask.h"
  30
  31 //__________________________________________________________________________________________________
  32 AliRsnAnalysisMonitorTask::AliRsnAnalysisMonitorTask(const char *name) :
  33   AliAnalysisTaskSE(name),
  34   fEventType(2),
  35   fNTracks(0),
  36   fOut(0x0),
  37   fTracks(0x0),
  38   fMaxITSband(1E6),
  39   fTPCpLimit(0.35),
  40   fLargeTPCband(-1E6),
  41   fSmallTPCband( 1E6),
  42   fESDtrackCutsTPC(),
  43   fESDtrackCutsITS(),
  44   fESDpid(0x0),
  45   fTOFmaker(0x0),
  46   fTOFcalib(0x0),
  47   fTOFcalibrateESD(kFALSE),
  48   fTOFcorrectTExp(kFALSE),
  49   fTOFuseT0(kFALSE),
  50   fTOFtuneMC(kFALSE),
  51   fTOFresolution(0.0)
  52
  53 {
  54 //
  55 // Constructor
  56 //
  57
  58   DefineOutput(1, TTree::Class());
  59 }
  60
  61 //__________________________________________________________________________________________________
  62 AliRsnAnalysisMonitorTask::AliRsnAnalysisMonitorTask(const AliRsnAnalysisMonitorTask& copy) :
  63   AliAnalysisTaskSE(copy),
  64   fEventType(2),
  65   fNTracks(0),
  66   fOut(0x0),
  67   fTracks(0x0),
  68   fMaxITSband(copy.fMaxITSband),
  69   fTPCpLimit(copy.fTPCpLimit),
  70   fLargeTPCband(copy.fLargeTPCband),
  71   fSmallTPCband(copy.fSmallTPCband),
  72   fESDtrackCutsTPC(copy.fESDtrackCutsTPC),
  73   fESDtrackCutsITS(copy.fESDtrackCutsITS),
  74   fESDpid(0x0),
  75   fTOFmaker(0x0),
  76   fTOFcalib(0x0),
  77   fTOFcalibrateESD(kFALSE),
  78   fTOFcorrectTExp(kFALSE),
  79   fTOFuseT0(kFALSE),
  80   fTOFtuneMC(kFALSE),
  81   fTOFresolution(0.0)
  82 {
  83 //
  84 // Copy constructor
  85 //
  86 }
  87
  88 //__________________________________________________________________________________________________
  89 AliRsnAnalysisMonitorTask& AliRsnAnalysisMonitorTask::operator=(const AliRsnAnalysisMonitorTask& copy)
  90 {
  91 //
  92 // Assignment operator
  93 //
  94
  95   fMaxITSband = copy.fMaxITSband;
  96
  97   fTPCpLimit    = copy.fTPCpLimit;
  98   fLargeTPCband = copy.fSmallTPCband;
  99   fSmallTPCband = copy.fLargeTPCband;
 100
 101   fESDtrackCutsTPC = copy.fESDtrackCutsTPC;
 102   fESDtrackCutsITS = copy.fESDtrackCutsITS;
 103
 104   fTOFcalibrateESD = copy.fTOFcalibrateESD;
 105   fTOFcorrectTExp = copy.fTOFcorrectTExp;
 106   fTOFuseT0 = copy.fTOFuseT0;
 107   fTOFtuneMC = copy.fTOFtuneMC;
 108   fTOFresolution = copy.fTOFresolution;
 109
 110   return (*this);
 111 }
 112
 113 //__________________________________________________________________________________________________
 114 AliRsnAnalysisMonitorTask::~AliRsnAnalysisMonitorTask()
 115 {
 116 //
 117 // Destructor
 118 //
 119
 120   if (fOut)    delete fOut;
 121   if (fESDpid) delete fESDpid;
 122 }
 123
 124 //__________________________________________________________________________________________________
 125 void AliRsnAnalysisMonitorTask::UserCreateOutputObjects()
 126 {
 127 //
 128 // Create the output data container
 129 //
 130
 131   // setup TPC response
 132   fESDpid = new AliESDpid;
 133   fESDpid->GetTPCResponse().SetBetheBlochParameters(fTPCpar[0], fTPCpar[1], fTPCpar[2], fTPCpar[3], fTPCpar[4]);
 134
 135   // setup TOF maker & calibration
 136   fTOFcalib = new AliTOFcalib;
 137   fTOFmaker = new AliTOFT0maker(fESDpid, fTOFcalib);
 138   fTOFmaker->SetTimeResolution(fTOFresolution);
 139
 140   // initialize all the branch arrays
 141   fTracks = new TClonesArray("AliRsnMonitorTrack", 0);
 142
 143   // create output tree
 144   OpenFile(1);
 145   fOut = new TTree("rsnMonitor", "Informations on single tracks for cut checking");
 146
 147   // link branches
 148   fOut->Branch("ntracks", &fNTracks     , "ntracks/I"   );
 149   fOut->Branch("evtype" , &fEventType   , "evtype/I"    );
 150   fOut->Branch("vertex" , &fVertex      , "vertex[3]/F");
 151   fOut->Branch("tracks" , "TClonesArray", &fTracks      );
 152 }
 153
 154 //__________________________________________________________________________________________________
 155 void AliRsnAnalysisMonitorTask::UserExec(Option_t *)
 156 {
 157 //
 158 // Main execution function.
 159 // Fills the fHEvents data member with the following legenda:
 160 // 0 -- event OK, prim vertex with tracks
 161 // 1 -- event OK, prim vertex with SPD
 162 // 2 -- event OK but vz large
 163 // 3 -- event bad
 164 //
 165
 166   static Int_t evNum = 0;
 167   evNum++;
 168
 169   // retrieve ESD event and related stack (if available)
 170   AliESDEvent *esd   = dynamic_cast<AliESDEvent*>(fInputEvent);
 171   AliStack    *stack = (fMCEvent ? fMCEvent->Stack() : 0x0);
 172
 173   // check the event
 174   EventEval(esd);
 175
 176   // if processable, then process it
 177   if      (fEventType == 0) ProcessESD(esd, esd->GetPrimaryVertexTracks(), stack);
 178   else if (fEventType == 1) ProcessESD(esd, esd->GetPrimaryVertexSPD()   , stack);
 179   else
 180   {
 181     fTracks->Delete();
 182     fTracks->Clear();
 183     fNTracks = 0;
 184   }
 185
 186   // add a new entry in the TTree
 187   fOut->Fill();
 188
 189   // update histogram container
 190   PostData(1, fOut);
 191 }
 192
 193 //__________________________________________________________________________________________________
 194 void AliRsnAnalysisMonitorTask::Terminate(Option_t *)
 195 {
 196 //
 197 // Terminate
 198 //
 199 }
 200
 201 //__________________________________________________________________________________________________
 202 void AliRsnAnalysisMonitorTask::EventEval(AliESDEvent *esd)
 203 {
 204 //
 205 // Checks if the event is good for analysis.
 206 // Sets the 'fEventType' flag to:
 207 // ---> 0 if a good primary vertex with tracks was found,
 208 // ---> 1 if a good SPD primary vertex was found
 209 // ---> 2 otherwise (event to be rejected)
 210 // In any case, adds an entry to the TTree, to keep trace of all events.
 211 //
 212
 213   // get number of tracks
 214   fNTracks = esd->GetNumberOfTracks();
 215
 216   // get the best primary vertex:
 217   // first try that with tracks, then the SPD one
 218   const AliESDVertex *vTrk  = esd->GetPrimaryVertexTracks();
 219   const AliESDVertex *vSPD  = esd->GetPrimaryVertexSPD();
 220   if(vTrk->GetNContributors() > 0)
 221   {
 222     fVertex[0] = vTrk->GetXv();
 223     fVertex[1] = vTrk->GetYv();
 224     fVertex[2] = vTrk->GetZv();
 225     fEventType = 0;
 226   }
 227   else if (vSPD->GetNContributors() > 0)
 228   {
 229     fVertex[0] = vSPD->GetXv();
 230     fVertex[1] = vSPD->GetYv();
 231     fVertex[2] = vSPD->GetZv();
 232     fEventType = 1;
 233   }
 234   else
 235   {
 236     fNTracks = 0;
 237     fEventType = 2;
 238   }
 239 }
 240
 241 //__________________________________________________________________________________________________
 242 void AliRsnAnalysisMonitorTask::ProcessESD
 243 (AliESDEvent *esd, const AliESDVertex *v, AliStack *stack)
 244 {
 245 //
 246 // Process the ESD container, to read all tracks and copy their useful values.
 247 // All info are stored into an AliRsnMonitorTrack object and saved into the
 248 // TClonesArray which is one of the branches of the output TTree.
 249 //
 250
 251   // clear array
 252   fTracks->Delete();
 253   fTracks->Clear();
 254
 255   // reject empty events
 256   if (!fNTracks) return;
 257
 258   // ITS stuff #1
 259   // create the response function and initialize it to MC or not
 260   // depending if the AliStack object is there or not
 261   Bool_t isMC = (stack != 0x0);
 262   AliITSPIDResponse itsrsp(isMC);
 263
 264   // TOF stuff #1: init OCDB
 265   Int_t run = esd->GetRunNumber();
 266   AliCDBManager *cdb = AliCDBManager::Instance();
 267   cdb->SetDefaultStorage("raw://");
 268   cdb->SetRun(run);
 269   // TOF stuff #2: init calibration
 270   fTOFcalib->SetCorrectTExp(fTOFcorrectTExp);
 271   fTOFcalib->Init();
 272   // TOF stuff #3: calibrate
 273   if (fTOFcalibrateESD) fTOFcalib->CalibrateESD(esd);
 274   if (fTOFtuneMC) fTOFmaker->TuneForMC(esd);
 275   if (fTOFuseT0)
 276   {
 277     fTOFmaker->ComputeT0TOF(esd);
 278     fTOFmaker->ApplyT0TOF(esd);
 279     fESDpid->MakePID(esd, kFALSE, 0.);
 280   }
 281   // TOF stuff #4: define fixed functions for compatibility range
 282   Double_t a1 = 0.01, a2 = -0.03;
 283   Double_t b1 = 0.25, b2 =  0.25;
 284   Double_t c1 = 0.05, c2 = -0.03;
 285   Double_t ymax, ymin;
 286
 287   // loop on all tracks
 288   Int_t               i, k, size, nITS;
 289   Double_t            tpcMaxNSigma, itsdedx[4], tofRef, tofRel;
 290   Bool_t              okTOF, okTrack, isTPC, isITSSA, matchedTOF;
 291   UChar_t             itsCluMap;
 292   Float_t             b[2], bCov[3];
 293   AliRsnMonitorTrack  mon;
 294
 295   for (i = 0; i < fNTracks; i++)
 296   {
 297     AliESDtrack *track = esd->GetTrack(i);
 298
 299     // skip NULL pointers, kink daughters and tracks which
 300     // cannot be propagated to primary vertex
 301     if (!track) continue;
 302     if ((Int_t)track->GetKinkIndex(0) > 0) continue;
 303     if (!track->RelateToVertex(v, esd->GetMagneticField(), kVeryBig)) continue;
 304
 305     // reset the output object
 306     // 'usable' flag will need to be set to 'ok'
 307     mon.Reset();
 308
 309     // get MC info if possible
 310     if (stack) mon.AdoptMC(TMath::Abs(track->GetLabel()), stack);
 311
 312     // copy general info
 313     mon.Status() = (UInt_t)track->GetStatus();
 314     mon.Length() = (Double_t)track->GetIntegratedLength();
 315     mon.Charge() = (Int_t)track->Charge();
 316     mon.PrecX()  = (Double_t)track->Px();
 317     mon.PrecY()  = (Double_t)track->Py();
 318     mon.PrecZ()  = (Double_t)track->Pz();
 319
 320     // evaluate some flags from the status to decide what to do next in some points
 321     isTPC      = ((mon.Status() & AliESDtrack::kTPCin)  != 0);
 322     isITSSA    = ((mon.Status() & AliESDtrack::kTPCin)  == 0 && (mon.Status() & AliESDtrack::kITSrefit) != 0 && (mon.Status() & AliESDtrack::kITSpureSA) == 0 && (mon.Status() & AliESDtrack::kITSpid) != 0);
 323     matchedTOF = ((mon.Status() & AliESDtrack::kTOFout) != 0 && (mon.Status() & AliESDtrack::kTIME) != 0);
 324
 325     // accept only tracks which are TPC+ITS or ITS standalone
 326     if (!isTPC && !isITSSA) continue;
 327
 328     // set the track type in the output object
 329     mon.ITSsa() = isITSSA;
 330
 331     // get DCA to primary vertex
 332     track->GetImpactParameters(b, bCov);
 333     mon.DCAr() = (Double_t)b[0];
 334     mon.DCAz() = (Double_t)b[1];
 335
 336     // get ITS info
 337     track->GetITSdEdxSamples(itsdedx);
 338     mon.ITSchi2() = track->GetITSchi2();
 339     mon.ITSsignal() = track->GetITSsignal();
 340     for (k = 0; k < 6; k++)
 341     {
 342       mon.ITSmap(k) = track->HasPointOnITSLayer(k);
 343       if (k < 4) mon.ITSdedx(k) = itsdedx[k];
 344     }
 345
 346     // get TPC info
 347     if (isTPC)
 348     {
 349       mon.TPCcount()  = (Int_t)track->GetTPCclusters(0);
 350       mon.TPCdedx()   = (Double_t)track->GetTPCsignal();
 351       mon.TPCchi2()   = (Double_t)track->GetTPCchi2();
 352       mon.TPCnsigma() = fESDpid->NumberOfSigmasTPC(track, AliPID::kKaon);
 353       mon.PtpcX()     = mon.PtpcY() = mon.PtpcZ() = 1E10;
 354       if (track->GetInnerParam())
 355       {
 356         mon.PtpcX() = track->GetInnerParam()->Px();
 357         mon.PtpcY() = track->GetInnerParam()->Py();
 358         mon.PtpcZ() = track->GetInnerParam()->Pz();
 359         for (k = 0; k < AliPID::kSPECIES; k++) mon.TPCref(k) = fESDpid->GetTPCResponse().GetExpectedSignal(mon.Ptpc(), (AliPID::EParticleType)k);
 360       }
 361     }
 362
 363     // get TOF info
 364     Double_t time[10];
 365     track->GetIntegratedTimes(time);
 366     mon.TOFsignal() = (Double_t)track->GetTOFsignal();
 367     for (k = 0; k < AliPID::kSPECIES; k++)
 368     {
 369       mon.TOFref(k)   = time[k];
 370       mon.TOFsigma(k) = (Double_t)fTOFmaker->GetExpectedSigma(mon.Prec(), time[k], AliPID::ParticleMass(k));
 371     }
 372
 373     // if we are here, the track is usable
 374     mon.SetUsable();
 375
 376     // now check the track against its cuts
 377     // and update the flag related to it
 378     // first, assume that cuts were passed
 379     // and if they aren't, just update the flag accordingly
 380     mon.CutsPassed() = kTRUE;
 381     if (TMath::Abs(mon.DCAz()) > 3.0) mon.CutsPassed() = kFALSE;
 382
 383     if (isTPC)
 384     {
 385       // check standard ESD cuts
 386       if (!fESDtrackCutsTPC.IsSelected(track)) mon.CutsPassed() = kFALSE;
 387
 388       // check TPC dE/dx
 389       if (mon.Ptpc() > fTPCpLimit) tpcMaxNSigma = fSmallTPCband; else tpcMaxNSigma = fLargeTPCband;
 390       if (TMath::Abs(mon.TPCnsigma()) > tpcMaxNSigma) mon.CutsPassed() = kFALSE;
 391
 392       // check TOF (only if momentum is large than function asymptote and flags are OK)
 393       okTOF = kTRUE;
 394       if (matchedTOF && mon.Prec() > TMath::Max(b1, b2))
 395       {
 396         tofRef = mon.TOFref(AliPID::kKaon);
 397         if (tofRef > 0.0)
 398         {
 399           tofRel   = (mon.TOFsignal() - tofRef) / tofRef;
 400           ymax     = a1 / (mon.Prec() - b1) + c1;
 401           ymin     = a2 / (mon.Prec() - b2) + c2;
 402           okTOF    = (tofRel >= ymin && tofRel <= ymax);
 403         }
 404       }
 405       if (!okTOF) mon.CutsPassed() = kFALSE;
 406     }
 407     else
 408     {
 409       // check standard ESD cuts
 410       if (!fESDtrackCutsITS.IsSelected(track)) mon.CutsPassed() = kFALSE;
 411
 412       // check dE/dx
 413       itsCluMap = track->GetITSClusterMap();
 414       nITS      = 0;
 415       for(k = 2; k < 6; k++) if(itsCluMap & (1 << k)) ++nITS;
 416       if (nITS < 3) okTrack = kFALSE; // track not good for PID
 417       mon.ITSnsigma() = itsrsp.GetNumberOfSigmas(mon.Prec(), mon.ITSsignal(), AliPID::kKaon, nITS, kTRUE);
 418       if (TMath::Abs(mon.ITSnsigma()) > fMaxITSband) mon.CutsPassed() = kFALSE;
 419     }
 420
 421     // collect only tracks which are declared usable
 422     if (mon.IsUsable())
 423     {
 424       size = (Int_t)fTracks->GetEntriesFast();
 425       new ((*fTracks)[size]) AliRsnMonitorTrack(mon);
 426     }
 427   }
 428 }