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->SetEtaMode(); //in case you want eta-Pt instead of y-Pt
  96   fProtonQAAnalysis->SetRunMCAnalysis();
  97   fProtonQAAnalysis->SetRunEfficiencyAnalysis(kFALSE); //kTRUE if the cuts should be used in the reco and pid efficiencies
  98   //fProtonQAAnalysis->SetMCProcessId(13);//4: weak decay - 13: hadronic interaction
  99   //fProtonQAAnalysis->SetMotherParticlePDGCode(3122);//3122: Lambda
 100
 101   //Use of TPConly tracks
 102   if(fProtonAnalysisMode == kTPC) {
 103     //fProtonQAAnalysis->SetQAYPtBins(20, -1.0, 1.0, 28, 0.1, 1.5); //TPC only
 104     fProtonQAAnalysis->SetQAYPtBins(10, -0.5, 0.5, 16, 0.5, 0.9); //TPC only
 105     fProtonQAAnalysis->UseTPCOnly();
 106     //fProtonQAAnalysis->SetTPCpid();
 107     fProtonQAAnalysis->SetMinTPCClusters(100);
 108     fProtonQAAnalysis->SetMaxChi2PerTPCCluster(2.2);
 109     fProtonQAAnalysis->SetMaxCov11(0.5);
 110     fProtonQAAnalysis->SetMaxCov22(0.5);
 111     fProtonQAAnalysis->SetMaxCov33(0.5);
 112     fProtonQAAnalysis->SetMaxCov44(0.5);
 113     fProtonQAAnalysis->SetMaxCov55(0.5);
 114     fProtonQAAnalysis->SetMaxSigmaToVertexTPC(2.0);
 115     //fProtonQAAnalysis->SetMaxDCAXYTPC(1.5);
 116     //fProtonQAAnalysis->SetMaxDCAZTPC(1.5);
 117   }
 118   //Use of HybridTPC tracks
 119   else if(fProtonAnalysisMode == kHybrid) {
 120     fProtonQAAnalysis->SetQAYPtBins(10, -0.5, 0.5, 16, 0.5, 0.9); //HybridTPC
 121     fProtonQAAnalysis->UseHybridTPC();
 122     fProtonQAAnalysis->SetTPCpid();
 123     fProtonQAAnalysis->SetMinTPCClusters(110);
 124     fProtonQAAnalysis->SetMaxChi2PerTPCCluster(2.2);
 125     fProtonQAAnalysis->SetMaxCov11(0.5);
 126     fProtonQAAnalysis->SetMaxCov22(0.5);
 127     fProtonQAAnalysis->SetMaxCov33(0.5);
 128     fProtonQAAnalysis->SetMaxCov44(0.5);
 129     fProtonQAAnalysis->SetMaxCov55(0.5);
 130     fProtonQAAnalysis->SetMaxSigmaToVertex(2.0);
 131     /*fProtonQAAnalysis->SetMaxDCAXY(1.5);
 132       fProtonQAAnalysis->SetMaxDCAZ(1.5);*/
 133     //fProtonQAAnalysis->SetPointOnITSLayer6();
 134     //fProtonQAAnalysis->SetPointOnITSLayer5();
 135     //fProtonQAAnalysis->SetPointOnITSLayer4();
 136     //fProtonQAAnalysis->SetPointOnITSLayer3();
 137     //fProtonQAAnalysis->SetPointOnITSLayer2();
 138     //fProtonQAAnalysis->SetPointOnITSLayer1();
 139     //fProtonQAAnalysis->SetMinITSClusters(4);
 140   }
 141   //Combined tracking
 142   else if(fProtonAnalysisMode == kGlobal) {
 143     fProtonQAAnalysis->SetQAYPtBins(20, -1.0, 1.0, 48, 0.3, 1.5); //combined tracking
 144     fProtonQAAnalysis->SetMinTPCClusters(110);
 145     fProtonQAAnalysis->SetMaxChi2PerTPCCluster(2.2);
 146     fProtonQAAnalysis->SetMaxCov11(0.5);
 147     fProtonQAAnalysis->SetMaxCov22(0.5);
 148     fProtonQAAnalysis->SetMaxCov33(0.5);
 149     fProtonQAAnalysis->SetMaxCov44(0.5);
 150     fProtonQAAnalysis->SetMaxCov55(0.5);
 151     fProtonQAAnalysis->SetMaxSigmaToVertex(2.0);
 152     //fProtonQAAnalysis->SetMaxDCAXY(2.0);
 153     //fProtonQAAnalysis->SetMaxDCAZ(2.0);
 154     fProtonQAAnalysis->SetTPCRefit();
 155     fProtonQAAnalysis->SetPointOnITSLayer1();
 156     fProtonQAAnalysis->SetPointOnITSLayer2();
 157     //fProtonQAAnalysis->SetPointOnITSLayer3();
 158     //fProtonQAAnalysis->SetPointOnITSLayer4();
 159     fProtonQAAnalysis->SetPointOnITSLayer5();
 160     fProtonQAAnalysis->SetPointOnITSLayer6();
 161     fProtonQAAnalysis->SetMinITSClusters(5);
 162     fProtonQAAnalysis->SetITSRefit();
 163     fProtonQAAnalysis->SetESDpid();
 164   }
 165
 166   fProtonQAAnalysis->SetPriorProbabilities(partFrac);
 167
 168   fList0 = new TList();
 169   fList0 = fProtonQAAnalysis->GetGlobalQAList();
 170
 171   fList1 = new TList();
 172   fList1 = fProtonQAAnalysis->GetPDGList();
 173
 174   fList2 = new TList();
 175   fList2 = fProtonQAAnalysis->GetMCProcessesList();
 176
 177   fList3 = new TList();
 178   fList3 = fProtonQAAnalysis->GetAcceptedCutList();
 179
 180   fList4 = new TList();
 181   fList4 = fProtonQAAnalysis->GetAcceptedDCAList();
 182
 183   fList5 = new TList();
 184   fList5 = fProtonQAAnalysis->GetEfficiencyQAList();
 185 }
 186
 187 //________________________________________________________________________
 188 void AliAnalysisTaskProtonsQA::Exec(Option_t *) {
 189   // Main loop
 190   // Called for each event
 191
 192   if (!fESD) {
 193     Printf("ERROR: fESD not available");
 194     return;
 195   }
 196
 197   if (!fMC) {
 198     Printf("ERROR: Could not retrieve MC event");
 199     return;
 200   }
 201
 202   AliStack* stack = fMC->Stack();
 203   if (!stack) {
 204     Printf("ERROR: Could not retrieve the stack");
 205     return;
 206   }
 207
 208   if(IsEventTriggered(fESD,fTriggerMode)) {
 209     const AliESDVertex *vertex = GetVertex(fESD,fProtonAnalysisMode,
 210                                            fVxMax,fVyMax,fVzMax);
 211     if(vertex) {
 212       fProtonQAAnalysis->RunQAAnalysis(stack, fESD, vertex);
 213       fProtonQAAnalysis->RunMCAnalysis(stack);
 214       fProtonQAAnalysis->RunEfficiencyAnalysis(stack, fESD, vertex);
 215     }//accepted vertex
 216   }//triggered event
 217
 218   // Post output data.
 219   PostData(0, fList0);
 220   PostData(1, fList1);
 221   PostData(2, fList2);
 222   PostData(3, fList3);
 223   PostData(4, fList4);
 224   PostData(5, fList5);
 225 }
 226
 227 //________________________________________________________________________
 228 void AliAnalysisTaskProtonsQA::Terminate(Option_t *) {
 229   // Draw result to the screen
 230   // Called once at the end of the query
 231
 232   fList0 = dynamic_cast<TList*> (GetOutputData(0));
 233   if (!fList0) {
 234     Printf("ERROR: fList0 not available");
 235     return;
 236   }
 237   fList1 = dynamic_cast<TList*> (GetOutputData(1));
 238   if (!fList1) {
 239     Printf("ERROR: fList1 not available");
 240     return;
 241   }
 242   fList2 = dynamic_cast<TList*> (GetOutputData(2));
 243   if (!fList2) {
 244     Printf("ERROR: fList2 not available");
 245     return;
 246   }
 247   fList3 = dynamic_cast<TList*> (GetOutputData(3));
 248   if (!fList3) {
 249     Printf("ERROR: fList3 not available");
 250     return;
 251   }
 252   fList4 = dynamic_cast<TList*> (GetOutputData(4));
 253   if (!fList4) {
 254     Printf("ERROR: fList4 not available");
 255     return;
 256   }
 257   fList5 = dynamic_cast<TList*> (GetOutputData(5));
 258   if (!fList5) {
 259     Printf("ERROR: fList5 not available");
 260     return;
 261   }
 262 }
 263
 264 //________________________________________________________________________
 265 Bool_t AliAnalysisTaskProtonsQA::IsEventTriggered(const AliESDEvent *esd,
 266                                                   TriggerMode trigger) {
 267   // check if the event was triggered
 268   ULong64_t triggerMask = esd->GetTriggerMask();
 269
 270   // definitions from p-p.cfg
 271   ULong64_t spdFO = (1 << 14);
 272   ULong64_t v0left = (1 << 11);
 273   ULong64_t v0right = (1 << 12);
 274
 275   switch (trigger) {
 276   case kMB1: {
 277     if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
 278       return kTRUE;
 279     break;
 280   }
 281   case kMB2: {
 282     if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
 283       return kTRUE;
 284     break;
 285   }
 286   case kSPDFASTOR: {
 287     if (triggerMask & spdFO)
 288       return kTRUE;
 289     break;
 290   }
 291   }//switch
 292
 293   return kFALSE;
 294 }
 295
 296 //________________________________________________________________________
 297 const AliESDVertex* AliAnalysisTaskProtonsQA::GetVertex(AliESDEvent* esd,
 298                                                         AnalysisMode mode,
 299                                                         Double_t gVxMax,
 300                                                         Double_t gVyMax,
 301                                                         Double_t gVzMax) {
 302   // Get the vertex from the ESD and returns it if the vertex is valid
 303   // Second argument decides which vertex is used (this selects
 304   // also the quality criteria that are applied)
 305   const AliESDVertex* vertex = 0;
 306   if(mode == kHybrid)
 307     vertex = esd->GetPrimaryVertexSPD();
 308   else if(mode == kTPC){
 309     Double_t kBz = esd->GetMagneticField();
 310     AliVertexerTracks vertexer(kBz);
 311     vertexer.SetTPCMode();
 312     AliESDVertex *vTPC = vertexer.FindPrimaryVertex(esd);
 313     esd->SetPrimaryVertexTPC(vTPC);
 314     for (Int_t i=0; i<esd->GetNumberOfTracks(); i++) {
 315       AliESDtrack *t = esd->GetTrack(i);
 316       t->RelateToVertexTPC(vTPC, kBz, kVeryBig);
 317     }
 318     delete vTPC;
 319     vertex = esd->GetPrimaryVertexTPC();
 320   }
 321   else if(mode == kGlobal)
 322     vertex = esd->GetPrimaryVertex();
 323   else
 324     Printf("GetVertex: ERROR: Invalid second argument %d", mode);
 325
 326   if(!vertex) return 0;
 327
 328   // check Ncontributors
 329   if(vertex->GetNContributors() <= 0) return 0;
 330
 331   // check resolution
 332   Double_t zRes = vertex->GetZRes();
 333   if(zRes == 0) return 0;
 334
 335   //check position
 336   if(TMath::Abs(vertex->GetXv()) > gVxMax) return 0;
 337   if(TMath::Abs(vertex->GetYv()) > gVyMax) return 0;
 338   if(TMath::Abs(vertex->GetZv()) > gVzMax) return 0;
 339
 340   return vertex;
 341 }
 342
 343
 344
 345