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