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 "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     }//accepted vertex
 237   }//triggered event
 238
 239   // Post output data.
 240   PostData(0, fList0);
 241   PostData(1, fList1);
 242   PostData(2, fList2);
 243   PostData(3, fList3);
 244   PostData(4, fList4);
 245   PostData(5, fList5);
 246   PostData(6, fList6);
 247   PostData(7, fList7);
 248 }
 249
 250 //________________________________________________________________________
 251 void AliAnalysisTaskProtonsQA::Terminate(Option_t *) {
 252   // Draw result to the screen
 253   // Called once at the end of the query
 254
 255   fList0 = dynamic_cast<TList*> (GetOutputData(0));
 256   if (!fList0) {
 257     Printf("ERROR: fList0 not available");
 258     return;
 259   }
 260   fList1 = dynamic_cast<TList*> (GetOutputData(1));
 261   if (!fList1) {
 262     Printf("ERROR: fList1 not available");
 263     return;
 264   }
 265   fList2 = dynamic_cast<TList*> (GetOutputData(2));
 266   if (!fList2) {
 267     Printf("ERROR: fList2 not available");
 268     return;
 269   }
 270   fList3 = dynamic_cast<TList*> (GetOutputData(3));
 271   if (!fList3) {
 272     Printf("ERROR: fList3 not available");
 273     return;
 274   }
 275   fList4 = dynamic_cast<TList*> (GetOutputData(4));
 276   if (!fList4) {
 277     Printf("ERROR: fList4 not available");
 278     return;
 279   }
 280   fList5 = dynamic_cast<TList*> (GetOutputData(5));
 281   if (!fList5) {
 282     Printf("ERROR: fList5 not available");
 283     return;
 284   }
 285 }
 286
 287 //________________________________________________________________________
 288 Bool_t AliAnalysisTaskProtonsQA::IsEventTriggered(const AliESDEvent *esd,
 289                                                   TriggerMode trigger) {
 290   // check if the event was triggered
 291   ULong64_t triggerMask = esd->GetTriggerMask();
 292
 293   // definitions from p-p.cfg
 294   ULong64_t spdFO = (1 << 14);
 295   ULong64_t v0left = (1 << 11);
 296   ULong64_t v0right = (1 << 12);
 297
 298   switch (trigger) {
 299   case kMB1: {
 300     if (triggerMask & spdFO || ((triggerMask & v0left) || (triggerMask & v0right)))
 301       return kTRUE;
 302     break;
 303   }
 304   case kMB2: {
 305     if (triggerMask & spdFO && ((triggerMask & v0left) || (triggerMask & v0right)))
 306       return kTRUE;
 307     break;
 308   }
 309   case kSPDFASTOR: {
 310     if (triggerMask & spdFO)
 311       return kTRUE;
 312     break;
 313   }
 314   }//switch
 315
 316   return kFALSE;
 317 }
 318
 319 //________________________________________________________________________
 320 const AliESDVertex* AliAnalysisTaskProtonsQA::GetVertex(AliESDEvent* esd,
 321                                                         AnalysisMode mode,
 322                                                         Double_t gVxMax,
 323                                                         Double_t gVyMax,
 324                                                         Double_t gVzMax) {
 325   // Get the vertex from the ESD and returns it if the vertex is valid
 326   // Second argument decides which vertex is used (this selects
 327   // also the quality criteria that are applied)
 328   const AliESDVertex* vertex = 0;
 329   if(mode == kHybrid)
 330     vertex = esd->GetPrimaryVertexSPD();
 331   else if(mode == kTPC){
 332     Double_t kBz = esd->GetMagneticField();
 333     AliVertexerTracks vertexer(kBz);
 334     vertexer.SetTPCMode();
 335     AliESDVertex *vTPC = vertexer.FindPrimaryVertex(esd);
 336     esd->SetPrimaryVertexTPC(vTPC);
 337     for (Int_t i=0; i<esd->GetNumberOfTracks(); i++) {
 338       AliESDtrack *t = esd->GetTrack(i);
 339       t->RelateToVertexTPC(vTPC, kBz, kVeryBig);
 340     }
 341     delete vTPC;
 342     vertex = esd->GetPrimaryVertexTPC();
 343   }
 344   else if(mode == kGlobal)
 345     vertex = esd->GetPrimaryVertex();
 346   else
 347     Printf("GetVertex: ERROR: Invalid second argument %d", mode);
 348
 349   if(!vertex) return 0;
 350
 351   // check Ncontributors
 352   if(vertex->GetNContributors() <= 0) return 0;
 353
 354   // check resolution
 355   Double_t zRes = vertex->GetZRes();
 356   if(zRes == 0) return 0;
 357
 358   //check position
 359   if(TMath::Abs(vertex->GetXv()) > gVxMax) return 0;
 360   if(TMath::Abs(vertex->GetYv()) > gVyMax) return 0;
 361   if(TMath::Abs(vertex->GetZv()) > gVzMax) return 0;
 362
 363   return vertex;
 364 }
 365
 366
 367
 368