PWGLF/RESONANCES/extra/AliAnalysisTaskResonanceQA.cxx

   1 #include <Riostream.h>
   2
   3 #include <TChain.h>
   4 #include <TH1.h>
   5 #include <TH2.h>
   6 #include <TH3.h>
   7 #include <TDatabasePDG.h>
   8 #include <TParticlePDG.h>
   9 #include <TParticle.h>
  10 #include <TLorentzVector.h>
  11
  12 #include "AliESDpid.h"
  13 #include "AliESDtrack.h"
  14 #include "AliESDEvent.h"
  15 #include "AliESDVertex.h"
  16 #include "AliESDtrackCuts.h"
  17
  18 #include "AliPID.h"
  19 #include "AliAnalysisTask.h"
  20 #include "AliAnalysisManager.h"
  21 #include "AliInputEventHandler.h"
  22
  23 #include "AliMCEvent.h"
  24 #include "AliStack.h"
  25 #include "AliMCParticle.h"
  26 #include "AliGenEventHeader.h"
  27
  28 #include "AliAnalysisTaskResonanceQA.h"
  29
  30 //
  31 // analysis task creating basic QA plots for resonance particles
  32 // Author: Ayben Karasu Uysal
  33 // Computes some QA histograms useful for checking productions and data
  34 // both for counting resonances inside MC, and for checking PID performances
  35 //
  36 // Revised by A. Pulvirenti
  37 //
  38
  39
  40 ClassImp(AliAnalysisTaskResonanceQA)
  41
  42 //________________________________________________________________________
  43 AliAnalysisTaskResonanceQA::AliAnalysisTaskResonanceQA(const char *name) :
  44    AliAnalysisTaskSE(name),
  45    fT0(AliESDpid::kTOF_T0),
  46    fPrimaryThr(1E-6),
  47    fVz(10.0),
  48    fOutputList(0),
  49    fSelectedEvts(0),
  50    fEventVz(0),
  51    fdEdxTPC(0),
  52    fdEdxITS(0),
  53    fTOFpid(0),
  54    fDCAXYvsPtBeforeCuts(0),
  55    fDCAZvsPtBeforeCuts(0),
  56    fNClusterPtBeforeCuts(0),
  57    fNFindableClusterPtBeforeCuts(0),
  58    fNCrossedRowsTPCPtBeforeCuts(0),
  59    fProducedParticles(0),
  60    fESD(0),
  61    fESDpid(0),
  62    fTrackCuts(0)
  63 {
  64 //
  65 // Constructor
  66 //
  67
  68    // Define input and output slots here
  69    // Input slot #0 works with a TChain
  70    DefineInput(0, TChain::Class());
  71
  72    // Output slot #0 id reserved by the base class for AOD
  73    // Output slot #1 writes into a TH1 container
  74    DefineOutput(1, TList::Class());
  75
  76    // setup to NULL all remaining histos
  77    Int_t i;
  78    for (i = 0; i < kResonances; i++) fRsnYPt[0][i] = fRsnYPt[1][i] = 0;
  79 }
  80
  81 //________________________________________________________________________
  82 const char* AliAnalysisTaskResonanceQA::EvtName(EEvtType type) const
  83 {
  84 //
  85 // Return a short string with name of resonance
  86 //
  87
  88    switch(type) {
  89       case kBadNoVertex   : return "NoVertex";
  90       case kBadVertexOutZ : return "FarVertex";
  91       case kGoodEvent     : return "Good";
  92       default             : return "Unknown";
  93    }
  94 }
  95
  96 //________________________________________________________________________
  97 const char* AliAnalysisTaskResonanceQA::RsnName(ERsn type) const
  98 {
  99 //
 100 // Return a short string with name of resonance
 101 //
 102
 103    switch(type) {
 104       case kPhi       : return "Phi";
 105       case kKStar0    : return "KStar0";
 106       case kRho       : return "Rho";
 107       case kLambdaStar: return "LambdaStar";
 108       case kXiStar0   : return "XiStar0";
 109       case kSigmaStarP: return "SigmaStarP";
 110       case kSigmaStarM: return "SigmaStarM";
 111       case kDeltaPP   : return "DeltaPP";
 112       default         : return "Unknown";
 113    }
 114 }
 115
 116 //________________________________________________________________________
 117 const char* AliAnalysisTaskResonanceQA::RsnSymbol(ERsn type) const
 118 {
 119 //
 120 // Return a short string with name of resonance
 121 //
 122
 123    switch(type) {
 124       case kPhi       : return "#phi";
 125       case kKStar0    : return "K^{*0}";
 126       case kRho       : return "#rho";
 127       case kLambdaStar: return "#Lambda^{*}";
 128       case kXiStar0   : return "#Xi^{*0}";
 129       case kSigmaStarP: return "#Sigma^{*+}";
 130       case kSigmaStarM: return "#Sigma^{*-}";
 131       case kDeltaPP   : return "#Delta^{++}";
 132       default         : return "Unknown";
 133    }
 134 }
 135
 136 //________________________________________________________________________
 137 Int_t AliAnalysisTaskResonanceQA::RsnPDG(ERsn type) const
 138 {
 139 //
 140 // Return PDG code of resonance
 141 //
 142
 143    switch(type) {
 144       case kPhi       : return  333;
 145       case kKStar0    : return  313;
 146       case kRho       : return  113;
 147       case kLambdaStar: return 3124;
 148       case kXiStar0   : return 3324;
 149       case kSigmaStarP: return 3224;
 150       case kSigmaStarM: return 3114;
 151       case kDeltaPP   : return 2224;
 152       default         : return    0;
 153    }
 154 }
 155
 156 //________________________________________________________________________
 157 void AliAnalysisTaskResonanceQA::UserCreateOutputObjects()
 158 {
 159 //
 160 // Create histograms (called once)
 161 //
 162
 163    Int_t i;
 164    Int_t vESDdEdxBin = 1000;  Double_t vESDdEdxMin = 0;  Double_t vESDdEdxMax = 1000;
 165    Int_t vESDQPBin = 900;     Double_t vESDQPMin = -4.5; Double_t vESDQPMax = 4.5;
 166    Int_t vPtBin = 500;        Double_t vPtMin = 0.;      Double_t vPtMax = 10.;
 167    Int_t vYBin = 600;         Double_t vYMin = -15.;     Double_t vYMax = 15.;
 168
 169    // standard objects
 170    fESDpid = new AliESDpid();
 171    fTrackCuts = AliESDtrackCuts::GetStandardITSTPCTrackCuts2010(kTRUE);
 172    fTrackCuts->SetPtRange(0.1, 1E20);
 173    fTrackCuts->SetEtaRange(-0.9, 0.9);
 174
 175    // create output list
 176    fOutputList = new TList();
 177    fOutputList->SetOwner();
 178
 179    // output histograms for PID signal
 180    fSelectedEvts = new TH1I("fSelectedEvts", "fSelectedEvts;fSelectedEvts", 4, 0, 4);
 181    fEventVz = new TH1F("fEventVz", "Vz distribution", 800, -40.0, 40.0);
 182    fdEdxTPC = new TH2F("fdEdxTPC", "dEdxTPC, After all cuts;chargexmomentum p (GeV/c); TPC signal (a.u.)", vESDQPBin, vESDQPMin, vESDQPMax, vESDdEdxBin, vESDdEdxMin, vESDdEdxMax);
 183    fdEdxITS = new TH2F("fdEdxITS", "dEdxITS, After all cuts;chargexmomentum p (GeV/c); TPC signal (a.u.)", vESDQPBin, vESDQPMin, vESDQPMax, vESDdEdxBin, vESDdEdxMin, vESDdEdxMax);
 184    fTOFpid = new TH2F("fTOFpid", "TOF PID;p (GeV/c);#beta;", vESDQPBin, vESDQPMin, vESDQPMax, 300, 0.1, 1.2);
 185
 186    // output histograms for track quality
 187    fDCAXYvsPtBeforeCuts = new TH2F("DCAXYvsPtBeforeCuts", "DCAXYvsPtBeforeCuts;p_{t} (GeV/c);DCAXY", vPtBin, vPtMin, vPtMax, 500, -10.0, 10.0);
 188    fDCAZvsPtBeforeCuts = new TH2F("DCAZvsPtBeforeCuts", "DCAZvsPtBeforeCuts;p_{t} (GeV/c);DCAZ", vPtBin, vPtMin, vPtMax, 500, -10.0, 10.0);
 189    fNClusterPtBeforeCuts = new TH2F("NClusterPtBeforeCuts", "NClusterPtBeforeCuts;p_{t} (GeV/c);NClusters", vPtBin, vPtMin, vPtMax, 180, 0., 180.);
 190    fNFindableClusterPtBeforeCuts = new TH2F("NFindableClusterPtBeforeCuts", "NFindableClusterPtBeforeCuts;p_{t} (GeV/c);NClusters", vPtBin, vPtMin, vPtMax, 180, 0., 180.);
 191    fNCrossedRowsTPCPtBeforeCuts = new TH2F("NCrossedRowsTPCPtBeforeCuts", "NCrossedRowsTPCPtBeforeCuts;p_{t} (GeV/c);NCrossedRowsTPC", vPtBin, vPtMin, vPtMax, 180, 0., 180.);
 192
 193    // output histograms for resonances
 194    fProducedParticles = new TH1I("ProducedParticles", "ProducedParticles;Produced Particles", kResonances, 0, kResonances);
 195    for (i = 0; i < kResonances; i++) {
 196       fRsnYPt[0][i] = new TH3F(Form("%s_all" , RsnName(i)), Form("%s -- ALL;p_{t} (GeV/c);Rapidity"    , RsnName(i)), vPtBin, vPtMin, vPtMax, vYBin, vYMin, vYMax, kEvtTypes, 0.0, (double)kEvtTypes);
 197       fRsnYPt[1][i] = new TH3F(Form("%s_prim", RsnName(i)), Form("%s -- PRIMARY;p_{t} (GeV/c);Rapidity", RsnName(i)), vPtBin, vPtMin, vPtMax, vYBin, vYMin, vYMax, kEvtTypes, 0.0, (double)kEvtTypes);
 198       fProducedParticles->GetXaxis()->SetBinLabel(i + 1, RsnSymbol(i));
 199    }
 200
 201    // add histogram to list
 202    fOutputList->Add(fSelectedEvts);
 203    fOutputList->Add(fEventVz);
 204    fOutputList->Add(fdEdxTPC);
 205    fOutputList->Add(fdEdxITS);
 206    fOutputList->Add(fTOFpid);
 207    fOutputList->Add(fDCAXYvsPtBeforeCuts);
 208    fOutputList->Add(fDCAZvsPtBeforeCuts);
 209    fOutputList->Add(fNClusterPtBeforeCuts);
 210    fOutputList->Add(fNFindableClusterPtBeforeCuts);
 211    fOutputList->Add(fNCrossedRowsTPCPtBeforeCuts);
 212    fOutputList->Add(fProducedParticles);
 213    for (i = 0; i < kResonances; i++) {
 214       fOutputList->Add(fRsnYPt[0][i]);
 215       fOutputList->Add(fRsnYPt[1][i]);
 216    }
 217
 218    PostData(1, fOutputList);
 219 }
 220
 221 //________________________________________________________________________
 222 void AliAnalysisTaskResonanceQA::UserExec(Option_t *)
 223 {
 224 //
 225 // Execute computations
 226 //
 227
 228    // first bin of this histogram counts processed events
 229    // second bin counts events passing physics selection
 230    // all events not passing physics selections are skipped
 231    fSelectedEvts->Fill(0.1);
 232    Bool_t isSelected = (((AliInputEventHandler*)(AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler()))->IsEventSelected() & AliVEvent::kMB);
 233    if (!isSelected) return;
 234
 235    // count CINT1B selected events
 236    fSelectedEvts->Fill(1.1);
 237
 238    // try to cast input event to ESD and loop on tracks
 239    fESD = dynamic_cast<AliESDEvent*>(InputEvent());
 240
 241    // check if event is accepted w.r. to all additional selection criteria
 242    Bool_t hasVertex = kFALSE, goodVertex = kFALSE;
 243    Double_t vz = 1E20, vzsign = 1E20;
 244    if (fESD) {
 245       // check primary vertex:
 246       // we accept only tracks/SPD with at least 1 contributor
 247       const AliESDVertex *v0 = fESD->GetPrimaryVertexTracks();
 248       if (v0) {
 249          if (v0->GetNContributors() > 0) {
 250             vzsign = v0->GetZ();
 251             vz = TMath::Abs(vzsign);
 252             hasVertex = kTRUE;
 253          } else {
 254             v0 = fESD->GetPrimaryVertexSPD();
 255             if (v0) {
 256                if (v0->GetNContributors() > 0) {
 257                   vzsign = v0->GetZ();
 258                   vz = TMath::Abs(vzsign);
 259                   hasVertex = kTRUE;
 260                }
 261             }
 262          }
 263       }
 264       if (hasVertex) {
 265          goodVertex = (vz <= fVz);
 266          fSelectedEvts->Fill(2.1);
 267          fEventVz->Fill(vzsign);
 268       }
 269    }
 270
 271    // determine event quality
 272    Bool_t eventAccepted = kFALSE;
 273    Int_t evtQuality = kEvtTypes;
 274    if (hasVertex) {
 275       if (goodVertex) {
 276          evtQuality = kGoodEvent;
 277          eventAccepted = kTRUE;
 278       } else {
 279          evtQuality = kBadVertexOutZ;
 280       }
 281    } else {
 282       evtQuality = kBadNoVertex;
 283    }
 284
 285    // Message
 286    static Int_t evNum = -1;
 287    evNum++;
 288    switch (evtQuality) {
 289       case kBadNoVertex  : AliInfo(Form("Rejecting event #%5d because it has not a good vertex", evNum)); break;
 290       case kBadVertexOutZ: AliInfo(Form("Rejecting event #%5d because it is out of range: vz = %5.3f vs. %5.3f", evNum, vz, fVz)); break;
 291       case kGoodEvent    : AliDebugClass(1, Form("Accepting event #%5d", evNum)); break;
 292       default            : AliError("This should never appear!!!");
 293    }
 294
 295    // if event is OK, loop on tracks to fill QA
 296    if (fESD && eventAccepted) {
 297       // if we arrive here, the event was processed successfully
 298       // then we fill last bin in first histogram
 299       fSelectedEvts->Fill(3.1);
 300
 301       // settings for TOF time zero
 302       if (fESD->GetTOFHeader())
 303          fESDpid->SetTOFResponse(fESD, AliESDpid::kTOF_T0);
 304       else {
 305          Float_t t0spread[10];
 306          Float_t intrinsicTOFres = 100;
 307          for (Int_t i = 0; i < 10; i++) t0spread[i] = (TMath::Sqrt(fESD->GetSigma2DiamondZ())) / 0.03; //0.03 to convert from cm to ps
 308          fESDpid->GetTOFResponse().SetT0resolution(t0spread);
 309          fESDpid->GetTOFResponse().SetTimeResolution(intrinsicTOFres);
 310          fESDpid->SetTOFResponse(fESD, fT0);
 311       }
 312       fESDpid->MakePID(fESD, kFALSE, 0.);
 313
 314       // use some variables to store useful info
 315       /*MISC   */ Int_t    iTrack, nTracks = fESD->GetNumberOfTracks();
 316       /*MISC   */ Double_t pt;
 317       /*QUALITY*/ Float_t  impactXY, impactZ, nCrossedRowsTPC;
 318       /*ITS    */ Double_t itsSignal = 0;
 319       /*TPC    */ Double_t ptotSE = 0, tpcSignal, signSE;
 320       /*TOF    */ Double_t momentum, length, time, beta; //, times[AliPID::kSPECIES], tzeroTrack = 0;
 321
 322       // loop on tracks
 323       for (iTrack = 0; iTrack < nTracks; iTrack++) {
 324
 325          // get track
 326          AliESDtrack* track = fESD->GetTrack(iTrack);
 327          if (!track) continue;
 328
 329          // get transverse momentum (used for quality histograms)
 330          pt = track->Pt();
 331
 332          // get charge and reject neutrals
 333          signSE = track->GetSign();
 334          if (signSE == 0) continue;
 335
 336          // get DCA and TPC-related quality params
 337          track->GetImpactParameters(impactXY, impactZ);
 338          nCrossedRowsTPC = track->GetTPCClusterInfo(2, 1);
 339
 340          // fill quality histograms when status flags are OK for TPC in+refit and ITS refit
 341          if (track->IsOn(AliESDtrack::kTPCin) && track->IsOn(AliESDtrack::kTPCrefit) && track->IsOn(AliESDtrack::kITSrefit)) {
 342             fNClusterPtBeforeCuts        ->Fill(pt, track->GetTPCNcls());
 343             fNFindableClusterPtBeforeCuts->Fill(pt, track->GetTPCNclsF());
 344             fNCrossedRowsTPCPtBeforeCuts ->Fill(pt, nCrossedRowsTPC);
 345             fDCAXYvsPtBeforeCuts         ->Fill(pt, impactXY);
 346             fDCAZvsPtBeforeCuts          ->Fill(pt, impactZ);
 347          }
 348
 349          // from now on, work only with tracks passing standard cuts for 2010
 350          if (!fTrackCuts->AcceptTrack(track)) continue;
 351
 352          // get momentum (used for PID histograms)
 353          momentum = track->P();
 354
 355          // ITS
 356          itsSignal = track->GetITSsignal();
 357          fdEdxITS->Fill(signSE * momentum, itsSignal);
 358
 359          // TPC
 360          if (track->GetInnerParam()) {
 361             AliExternalTrackParam trackIn1(*track->GetInnerParam());
 362             ptotSE = trackIn1.GetP();
 363             tpcSignal = track->GetTPCsignal();
 364             fdEdxTPC->Fill(signSE * ptotSE, tpcSignal);
 365          } // end TPC
 366
 367          // TOF (requires checking some more status flags
 368          if (track->IsOn(AliESDtrack::kTOFpid) && track->IsOn(AliESDtrack::kTOFout) && track->IsOn(AliESDtrack::kTIME) && !track->IsOn(AliESDtrack::kTOFmismatch)) {
 369             length = track->GetIntegratedLength();
 370
 371             // reject suspiciously small lengths/times
 372             if (track->GetIntegratedLength() < 350.) continue;
 373             if (track->GetTOFsignal() < 1E-6) continue;
 374
 375             // thhese are maybe not needed
 376             // track->GetIntegratedTimes(times);
 377             // tzeroTrack = fESDpid->GetTOFResponse().GetStartTime(track->P());
 378
 379             // get TOF measured time and compute beta
 380             time = track->GetTOFsignal() - fESDpid->GetTOFResponse().GetStartTime(momentum);
 381             beta = length / (2.99792457999999984e-02 * time);
 382             fTOFpid->Fill(signSE * momentum, beta);
 383          } // end TOF
 384       } // end track loop
 385    } // end ESD check
 386
 387    // loop on MC, if MC event and stack are available
 388    AliMCEvent *mcEvent = MCEvent();
 389    if (mcEvent) {
 390       // MC primary vertex
 391       TArrayF mcVertex(3);
 392       AliGenEventHeader *genEH = mcEvent->GenEventHeader();
 393       genEH->PrimaryVertex(mcVertex);
 394       // loop on particles
 395       AliStack *stack = mcEvent->Stack();
 396       if (stack) {
 397          Int_t iTrack, irsn, pdg, mother, motherPDG, nTracks = stack->GetNtrack();
 398          Double_t dx, dy, dz, dist;
 399          TLorentzVector vprod;
 400          for (iTrack = 0; iTrack < nTracks; iTrack++) {
 401
 402             // get particle
 403             TParticle *part = stack->Particle(iTrack);
 404             if (!part) {
 405                AliError(Form("Could not receive track %d", iTrack));
 406                continue;
 407             }
 408
 409             // get PDG code of particle and check physical primary
 410             pdg = part->GetPdgCode();
 411             mother = part->GetFirstMother();
 412             motherPDG = 0;
 413
 414             // loop over possible resonances
 415             for (irsn = 0; irsn < kResonances; irsn++) {
 416                if (TMath::Abs(pdg) == RsnPDG(irsn)) {
 417                   // debug check
 418                   if (mother >= 0) {
 419                      TParticle *p = stack->Particle(mother);
 420                      motherPDG = p->GetPdgCode();
 421                   }
 422                   part->ProductionVertex(vprod);
 423                   dx = mcVertex[0] - vprod.X();
 424                   dy = mcVertex[1] - vprod.Y();
 425                   dz = mcVertex[2] - vprod.Z();
 426                   dist = TMath::Sqrt(dx*dx + dy*dy + dz*dz);
 427                   AliDebugClass(1, Form("PDG = %d -- mother ID = %d, pdg = %d -- dist. from MC vertex = %f -- prod. time = %f", pdg, mother, motherPDG, dist, vprod.T()));
 428                   // global counter is always filled
 429                   fRsnYPt[0][irsn]->Fill(part->Pt(), part->Y(), ((double)(evtQuality) + 0.2));
 430                   // counter for primaries is filled only if mother is less than 0 (primary)
 431                   if (dist < fPrimaryThr) {
 432                      fRsnYPt[1][irsn]->Fill(part->Pt(), part->Y(), ((double)(evtQuality) + 0.2));
 433                   }
 434                   // fill the global histogram
 435                   fProducedParticles->Fill(irsn + 1);
 436                   // if one PDG match was found, no need to continue loop
 437                   break;
 438                }
 439             }
 440          }
 441       }
 442    }
 443
 444    PostData(1, fOutputList);
 445 }
 446
 447 //________________________________________________________________________
 448 void AliAnalysisTaskResonanceQA::Terminate(Option_t *)
 449 {
 450 //
 451 // Termination
 452 // Quickly check that the output list is there
 453 //
 454
 455    fOutputList = dynamic_cast<TList*>(GetOutputData(1));
 456    if (!fOutputList) {
 457       AliError("Output list not found!!!");
 458       return;
 459    }
 460 }