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