PWG2/AliAnalysisTaskProtonsQA.cxx

   1 #include "TChain.h"
   2 #include "TTree.h"
   3 #include "TString.h"
   4 #include "TList.h"
   5 #include "TH2F.h"
   6 #include "TH1I.h"
   7 #include "TF1.h"
   8 #include "TCanvas.h"
   9 #include "TLorentzVector.h"
  10
  11 #include "AliAnalysisTask.h"
  12 #include "AliAnalysisManager.h"
  13
  14 #include "AliESDEvent.h"
  15 #include "AliESDInputHandler.h"
  16 #include "AliMCEventHandler.h"
  17 #include "AliMCEvent.h"
  18 #include "AliStack.h"
  19 #include "AliVertexerTracks.h"
  20 #include "AliESDVertex.h"
  21
  22 #include "AliProtonQAAnalysis.h"
  23 #include "AliAnalysisTaskProtonsQA.h"
  24
  25 // Analysis task used for the QA of the (anti)proton analysis
  26 // Author: Panos Cristakoglou
  27
  28 ClassImp(AliAnalysisTaskProtonsQA)
  29
  30 //________________________________________________________________________
  31 AliAnalysisTaskProtonsQA::AliAnalysisTaskProtonsQA()
  32   : AliAnalysisTask(), fESD(0), fMC(0),
  33     fList0(0), fList1(0), fList2(0), fList3(0), fList4(0), fList5(0),
  34     fProtonQAAnalysis(0),
  35     fTriggerMode(kMB2), fProtonAnalysisMode(kTPC),
  36     fVxMax(0), fVyMax(0), fVzMax(0) {
  37   //Dummy constructor
  38 }
  39
  40 //________________________________________________________________________
  41 AliAnalysisTaskProtonsQA::AliAnalysisTaskProtonsQA(const char *name)
  42   : AliAnalysisTask(name, ""), fESD(0), fMC(0),
  43     fList0(0), fList1(0), fList2(0), fList3(0), fList4(0), fList5(0),
  44     fProtonQAAnalysis(0),
  45     fTriggerMode(kMB2), fProtonAnalysisMode(kTPC),
  46     fVxMax(0), fVyMax(0), fVzMax(0) {
  47   // Constructor
  48
  49   // Define input and output slots here
  50   // Input slot #0 works with a TChain
  51   DefineInput(0, TChain::Class());
  52   // Output slot #0 writes into a TList container
  53   DefineOutput(0, TList::Class());
  54   DefineOutput(1, TList::Class());
  55   DefineOutput(2, TList::Class());
  56   DefineOutput(3, TList::Class());
  57   DefineOutput(4, TList::Class());
  58   DefineOutput(5, TList::Class());
  59 }
  60
  61 //________________________________________________________________________
  62 void AliAnalysisTaskProtonsQA::ConnectInputData(Option_t *) {
  63   // Connect ESD or AOD here
  64   // Called once
  65
  66   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
  67   if (!tree) {
  68     Printf("ERROR: Could not read chain from input slot 0");
  69   } else {
  70     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
  71
  72     if (!esdH) {
  73       Printf("ERROR: Could not get ESDInputHandler");
  74     } else
  75       fESD = esdH->GetEvent();
  76   }
  77
  78   AliMCEventHandler* mcH = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
  79   if (!mcH) {
  80     Printf("ERROR: Could not retrieve MC event handler");
  81   }
  82   else
  83     fMC = mcH->MCEvent();
  84 }
  85
  86 //________________________________________________________________________
  87 void AliAnalysisTaskProtonsQA::CreateOutputObjects() {
  88   // Create histograms
  89   // Called once
  90   //Prior probabilities
  91   Double_t partFrac[5] = {0.01, 0.01, 0.85, 0.10, 0.05};
  92
  93   //proton analysis object
  94   fProtonQAAnalysis = new AliProtonQAAnalysis();
  95   fProtonQAAnalysis->SetRunMCAnalysis();
  96   fProtonQAAnalysis->SetRunEfficiencyAnalysis(kTRUE,kFALSE); //kTRUE,kTRUE for eta-pT efficiencies and if the cuts should be used in the reco and pid efficiencies
  97   //fProtonQAAnalysis->SetMCProcessId(13);//4: weak decay - 13: hadronic interaction
  98   //fProtonQAAnalysis->SetMotherParticlePDGCode(3122);//3122: Lambda
  99
 100   //Use of TPConly tracks
 101   if(fProtonAnalysisMode == kTPC) {
 102     fProtonQAAnalysis->SetQAYPtBins(10, -0.5, 0.5, 12, 0.5, 0.9); //TPC only
 103     fProtonQAAnalysis->UseTPCOnly();
 104     //fProtonQAAnalysis->SetTPCpid();
 105     fProtonQAAnalysis->SetMinTPCClusters(100);
 106     fProtonQAAnalysis->SetMaxChi2PerTPCCluster(2.2);
 107     fProtonQAAnalysis->SetMaxCov11(0.5);
 108     fProtonQAAnalysis->SetMaxCov22(0.5);
 109     fProtonQAAnalysis->SetMaxCov33(0.5);
 110     fProtonQAAnalysis->SetMaxCov44(0.5);
 111     fProtonQAAnalysis->SetMaxCov55(0.5);
 112     fProtonQAAnalysis->SetMaxSigmaToVertexTPC(2.0);
 113     //fProtonQAAnalysis->SetMaxDCAXYTPC(1.5);
 114     //fProtonQAAnalysis->SetMaxDCAZTPC(1.5);
 115   }
 116   //Use of HybridTPC tracks
 117   else if(fProtonAnalysisMode == kHybrid) {
 118     fProtonQAAnalysis->SetQAYPtBins(10, -0.5, 0.5, 12, 0.5, 0.9); //HybridTPC
 119     fProtonQAAnalysis->UseHybridTPC();
 120     fProtonQAAnalysis->SetTPCpid();
 121     fProtonQAAnalysis->SetMinTPCClusters(110);
 122     fProtonQAAnalysis->SetMaxChi2PerTPCCluster(2.2);
 123     fProtonQAAnalysis->SetMaxCov11(0.5);
 124     fProtonQAAnalysis->SetMaxCov22(0.5);
 125     fProtonQAAnalysis->SetMaxCov33(0.5);
 126     fProtonQAAnalysis->SetMaxCov44(0.5);
 127     fProtonQAAnalysis->SetMaxCov55(0.5);
 128     fProtonQAAnalysis->SetMaxSigmaToVertex(2.0);
 129     /*fProtonQAAnalysis->SetMaxDCAXY(1.5);
 130       fProtonQAAnalysis->SetMaxDCAZ(1.5);*/
 131     fProtonQAAnalysis->SetPointOnITSLayer6();
 132     fProtonQAAnalysis->SetPointOnITSLayer5();
 133   //fProtonQAAnalysis->SetPointOnITSLayer4();
 134   //fProtonQAAnalysis->SetPointOnITSLayer3();
 135     fProtonQAAnalysis->SetPointOnITSLayer2();
 136     fProtonQAAnalysis->SetPointOnITSLayer1();
 137     fProtonQAAnalysis->SetMinITSClusters(5);
 138   }
 139   //Combined tracking
 140   else if(fProtonAnalysisMode == kGlobal) {
 141     fProtonQAAnalysis->SetQAYPtBins(20, -1.0, 1.0, 27, 0.4, 3.1); //combined tracking
 142     fProtonQAAnalysis->SetMinTPCClusters(110);
 143     fProtonQAAnalysis->SetMaxChi2PerTPCCluster(2.2);
 144     fProtonQAAnalysis->SetMaxCov11(0.5);
 145     fProtonQAAnalysis->SetMaxCov22(0.5);
 146     fProtonQAAnalysis->SetMaxCov33(0.5);
 147     fProtonQAAnalysis->SetMaxCov44(0.5);
 148     fProtonQAAnalysis->SetMaxCov55(0.5);
 149     fProtonQAAnalysis->SetMaxSigmaToVertex(2.0);
 150     //fProtonQAAnalysis->SetMaxDCAXY(2.0);
 151     //fProtonQAAnalysis->SetMaxDCAZ(2.0);
 152     fProtonQAAnalysis->SetTPCRefit();
 153     fProtonQAAnalysis->SetPointOnITSLayer1();
 154     fProtonQAAnalysis->SetPointOnITSLayer2();
 155     //fProtonQAAnalysis->SetPointOnITSLayer3();
 156     //fProtonQAAnalysis->SetPointOnITSLayer4();
 157     fProtonQAAnalysis->SetPointOnITSLayer5();
 158     fProtonQAAnalysis->SetPointOnITSLayer6();
 159     fProtonQAAnalysis->SetMinITSClusters(5);
 160     fProtonQAAnalysis->SetITSRefit();
 161     fProtonQAAnalysis->SetESDpid();
 162   }
 163
 164   fProtonQAAnalysis->SetPriorProbabilities(partFrac);
 165
 166   fList0 = new TList();
 167   fList0 = fProtonQAAnalysis->GetGlobalQAList();
 168
 169   fList1 = new TList();
 170   fList1 = fProtonQAAnalysis->GetPDGList();
 171
 172   fList2 = new TList();
 173   fList2 = fProtonQAAnalysis->GetMCProcessesList();
 174
 175   fList3 = new TList();
 176   fList3 = fProtonQAAnalysis->GetAcceptedCutList();
 177
 178   fList4 = new TList();
 179   fList4 = fProtonQAAnalysis->GetAcceptedDCAList();
 180
 181   fList5 = new TList();
 182   fList5 = fProtonQAAnalysis->GetEfficiencyQAList();
 183 }
 184
 185 //________________________________________________________________________
 186 void AliAnalysisTaskProtonsQA::Exec(Option_t *) {
 187   // Main loop
 188   // Called for each event
 189
 190   if (!fESD) {
 191     Printf("ERROR: fESD not available");
 192     return;
 193   }
 194
 195   if (!fMC) {
 196     Printf("ERROR: Could not retrieve MC event");
 197     return;
 198   }
 199
 200   AliStack* stack = fMC->Stack();
 201   if (!stack) {
 202     Printf("ERROR: Could not retrieve the stack");
 203     return;
 204   }
 205
 206   if(IsEventTriggered(fESD,fTriggerMode)) {
 207     const AliESDVertex *vertex = GetVertex(fESD,fProtonAnalysisMode,
 208                                            fVxMax,fVyMax,fVzMax);
 209     if(vertex) {
 210       fProtonQAAnalysis->RunQAAnalysis(stack, fESD);
 211       fProtonQAAnalysis->RunMCAnalysis(stack);
 212       fProtonQAAnalysis->RunEfficiencyAnalysis(stack, fESD);
 213     }//accepted vertex
 214   }//triggered event
 215
 216   // Post output data.
 217   PostData(0, fList0);
 218   PostData(1, fList1);
 219   PostData(2, fList2);
 220   PostData(3, fList3);
 221   PostData(4, fList4);
 222   PostData(5, fList5);
 223 }
 224
 225 //________________________________________________________________________
 226 void AliAnalysisTaskProtonsQA::Terminate(Option_t *) {
 227   // Draw result to the screen
 228   // Called once at the end of the query
 229
 230   fList0 = dynamic_cast<TList*> (GetOutputData(0));
 231   if (!fList0) {
 232     Printf("ERROR: fList0 not available");
 233     return;
 234   }
 235   fList1 = dynamic_cast<TList*> (GetOutputData(1));
 236   if (!fList1) {
 237     Printf("ERROR: fList1 not available");
 238     return;
 239   }
 240   fList2 = dynamic_cast<TList*> (GetOutputData(2));
 241   if (!fList2) {
 242     Printf("ERROR: fList2 not available");
 243     return;
 244   }
 245   fList3 = dynamic_cast<TList*> (GetOutputData(3));
 246   if (!fList3) {
 247     Printf("ERROR: fList3 not available");
 248     return;
 249   }
 250   fList4 = dynamic_cast<TList*> (GetOutputData(4));
 251   if (!fList4) {
 252     Printf("ERROR: fList4 not available");
 253     return;
 254   }
 255   fList5 = dynamic_cast<TList*> (GetOutputData(5));
 256   if (!fList5) {
 257     Printf("ERROR: fList5 not available");
 258     return;
 259   }
 260 }
 261
 262 //________________________________________________________________________
 263 Bool_t AliAnalysisTaskProtonsQA::IsEventTriggered(const AliESDEvent *esd,
 264                                                   TriggerMode trigger) {
 265   // check if the event was triggered
 266   ULong64_t triggerMask = esd->GetTriggerMask();
 267
 268   // definitions from p-p.cfg
 269   ULong64_t spdFO = (1 << 14);
 270   ULong64_t v0left = (1 << 11);
 271   ULong64_t v0right = (1 << 12);
 272
 273   switch (trigger) {
 274   case kMB1: {
 275     if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
 276       return kTRUE;
 277     break;
 278   }
 279   case kMB2: {
 280     if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
 281       return kTRUE;
 282     break;
 283   }
 284   case kSPDFASTOR: {
 285     if (triggerMask & spdFO)
 286       return kTRUE;
 287     break;
 288   }
 289   }//switch
 290
 291   return kFALSE;
 292 }
 293
 294 //________________________________________________________________________
 295 const AliESDVertex* AliAnalysisTaskProtonsQA::GetVertex(AliESDEvent* esd,
 296                                                         AnalysisMode mode,
 297                                                         Double_t gVxMax,
 298                                                         Double_t gVyMax,
 299                                                         Double_t gVzMax) {
 300   // Get the vertex from the ESD and returns it if the vertex is valid
 301   // Second argument decides which vertex is used (this selects
 302   // also the quality criteria that are applied)
 303   const AliESDVertex* vertex = 0;
 304   if(mode == kHybrid)
 305     vertex = esd->GetPrimaryVertexSPD();
 306   else if(mode == kTPC){
 307     Double_t kBz = esd->GetMagneticField();
 308     AliVertexerTracks vertexer(kBz);
 309     vertexer.SetTPCMode();
 310     AliESDVertex *vTPC = vertexer.FindPrimaryVertex(esd);
 311     esd->SetPrimaryVertexTPC(vTPC);
 312     for (Int_t i=0; i<esd->GetNumberOfTracks(); i++) {
 313       AliESDtrack *t = esd->GetTrack(i);
 314       t->RelateToVertexTPC(vTPC, kBz, kVeryBig);
 315     }
 316     delete vTPC;
 317     vertex = esd->GetPrimaryVertexTPC();
 318   }
 319   else if(mode == kGlobal)
 320     vertex = esd->GetPrimaryVertex();
 321   else
 322     Printf("GetVertex: ERROR: Invalid second argument %d", mode);
 323
 324   if(!vertex) return 0;
 325
 326   // check Ncontributors
 327   if(vertex->GetNContributors() <= 0) return 0;
 328
 329   // check resolution
 330   Double_t zRes = vertex->GetZRes();
 331   if(zRes == 0) return 0;
 332
 333   //check position
 334   if(TMath::Abs(vertex->GetXv()) > gVxMax) return 0;
 335   if(TMath::Abs(vertex->GetYv()) > gVyMax) return 0;
 336   if(TMath::Abs(vertex->GetZv()) > gVzMax) return 0;
 337
 338   return vertex;
 339 }
 340
 341
 342
 343