HMPID/AliHMPIDAnalysisTask.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //==============================================================================
  17 // AliHMPIDAnalysysTask - Class representing a basic analysis tool of HMPID data at
  18 // level of ESD.
  19 // A set of histograms is created.
  20 //==============================================================================
  21 //
  22 // By means of AliHMPIDAnalysisTask.C macro it is possible to use this class
  23 // to perform the analysis on local data, on data on alien using local machine
  24 // and on CAF.
  25
  26 #ifndef AliHMPIDAnalysisTASK_CXX
  27 #define AliHMPIDAnalysisTASK_CXX
  28
  29
  30 #include "TH1.h"
  31 #include "TH2.h"
  32 #include "TFile.h"
  33 #include "TCanvas.h"
  34 #include "TGraphErrors.h"
  35 #include "AliAnalysisManager.h"
  36 #include "AliESDInputHandler.h"
  37 #include "AliMCEventHandler.h"
  38 #include "AliMCEvent.h"
  39 #include "AliESDtrack.h"
  40 #include "AliPID.h"
  41 #include "AliLog.h"
  42 #include "AliHMPIDAnalysisTask.h"
  43
  44 ClassImp(AliHMPIDAnalysisTask)
  45
  46 //__________________________________________________________________________
  47 AliHMPIDAnalysisTask::AliHMPIDAnalysisTask() :
  48   fESD(0x0),fHmpHistList(0x0),
  49   fNevts(0),
  50   fTrigNevts(0),
  51   fTrigger(0),
  52   fHmpPesdPhmp(0x0),fHmpCkovPesd(0x0),fHmpCkovPhmp(0x0),
  53   fHmpMipTrkDist(0x0),fHmpMipTrkDistX(0x0),fHmpMipTrkDistY(0x0),fHmpMipCharge3cm(0x0),fHmpMipCharge1cm(0x0),fHmpNumPhots(0x0),
  54   fHmpTrkFlags(0x0),
  55   fN1(6),
  56   fN2(8),
  57   fPionEff(0x0),
  58   fKaonEff(0x0),
  59   fProtEff(0x0),
  60   fPionTot(0x0),
  61   fKaonTot(0x0),
  62   fProtTot(0x0),
  63   fPionNot(0x0),
  64   fKaonNot(0x0),
  65   fProtNot(0x0),
  66   fPionCon(0x0),
  67   fKaonCon(0x0),
  68   fProtCon(0x0),
  69   fThetavsPiFromK(0x0),
  70   fThetapivsPesd(0x0),
  71   fThetaKvsPesd(0x0),
  72   fThetaPvsPesd(0x0)
  73 {
  74   //
  75   //Default ctor
  76   //
  77 }
  78
  79 //___________________________________________________________________________
  80 AliHMPIDAnalysisTask::AliHMPIDAnalysisTask(const Char_t* name) :
  81   AliAnalysisTaskSE(name),
  82   fESD(0), fHmpHistList(0x0), fNevts(0),
  83   fTrigNevts(0),
  84   fTrigger(0),
  85   fHmpPesdPhmp(0x0), fHmpCkovPesd(0x0), fHmpCkovPhmp(0x0),
  86   fHmpMipTrkDist(0x0), fHmpMipTrkDistX(0x0), fHmpMipTrkDistY(0x0),
  87   fHmpMipCharge3cm(0x0), fHmpMipCharge1cm(0x0),
  88   fHmpNumPhots(0x0), fHmpTrkFlags(0x0),
  89   fN1(6),
  90   fN2(8),
  91   fPionEff(0x0),
  92   fKaonEff(0x0),
  93   fProtEff(0x0),
  94   fPionTot(0x0),
  95   fKaonTot(0x0),
  96   fProtTot(0x0),
  97   fPionNot(0x0),
  98   fKaonNot(0x0),
  99   fProtNot(0x0),
 100   fPionCon(0x0),
 101   fKaonCon(0x0),
 102   fProtCon(0x0),
 103   fThetavsPiFromK(0x0),
 104   fThetapivsPesd(0x0),
 105   fThetaKvsPesd(0x0),
 106   fThetaPvsPesd(0x0)
 107 {
 108   //
 109   // Constructor. Initialization of Inputs and Outputs
 110   //
 111
 112   DefineOutput(0,TList::Class());
 113 }
 114
 115 //___________________________________________________________________________
 116 AliHMPIDAnalysisTask& AliHMPIDAnalysisTask::operator=(const AliHMPIDAnalysisTask& c)
 117 {
 118   //
 119   // Assignment operator
 120   //
 121   if (this!=&c) {
 122     AliAnalysisTaskSE::operator=(c);
 123     fESD             = c.fESD;
 124     fHmpHistList     = c.fHmpHistList;
 125     fNevts           = c.fNevts;
 126     fTrigNevts       = c.fTrigNevts;
 127     fTrigger         = c.fTrigger;
 128     fHmpPesdPhmp     = c.fHmpPesdPhmp;
 129     fHmpCkovPesd     = c.fHmpCkovPesd;
 130     fHmpCkovPhmp     = c.fHmpCkovPhmp;
 131     fHmpMipTrkDist   = c.fHmpMipTrkDist;
 132     fHmpMipTrkDistX  = c.fHmpMipTrkDistX;
 133     fHmpMipTrkDistY  = c.fHmpMipTrkDistY;
 134     fHmpMipCharge3cm = c.fHmpMipCharge3cm;
 135     fHmpMipCharge1cm = c.fHmpMipCharge1cm;
 136     fHmpNumPhots     = c.fHmpNumPhots;
 137     fHmpTrkFlags     = c.fHmpTrkFlags;
 138     fN1              = c.fN1;
 139     fN2              = c.fN2;
 140     fPionEff         = c.fPionEff;
 141     fKaonEff         = c.fKaonEff;
 142     fProtEff         = c.fProtEff;
 143     fPionTot         = c.fPionTot;
 144     fKaonTot         = c.fKaonTot;
 145     fProtTot         = c.fProtTot;
 146     fPionNot         = c.fPionNot;
 147     fKaonNot         = c.fKaonNot;
 148     fProtNot         = c.fProtNot;
 149     fPionCon         = c.fPionCon;
 150     fKaonCon         = c.fKaonCon;
 151     fProtCon         = c.fProtCon;
 152     fThetavsPiFromK  = c.fThetavsPiFromK;
 153     fThetapivsPesd   = c.fThetapivsPesd;
 154     fThetaKvsPesd    = c.fThetaKvsPesd;
 155     fThetaPvsPesd    = c.fThetaPvsPesd;
 156   }
 157   return *this;
 158 }
 159
 160 //___________________________________________________________________________
 161 AliHMPIDAnalysisTask::AliHMPIDAnalysisTask(const AliHMPIDAnalysisTask& c) :
 162   AliAnalysisTaskSE(c),
 163   fESD(c.fESD),fHmpHistList(c.fHmpHistList), fNevts(c.fNevts),
 164   fTrigNevts(c.fTrigNevts),
 165   fTrigger(c.fTrigger),
 166   fHmpPesdPhmp(c.fHmpPesdPhmp),fHmpCkovPesd(c.fHmpCkovPesd),fHmpCkovPhmp(c.fHmpCkovPhmp),
 167   fHmpMipTrkDist(c.fHmpMipTrkDist),fHmpMipTrkDistX(c.fHmpMipTrkDistX),fHmpMipTrkDistY(c.fHmpMipTrkDistY),
 168   fHmpMipCharge3cm(c.fHmpMipCharge3cm), fHmpMipCharge1cm(c.fHmpMipCharge1cm),
 169   fHmpNumPhots(c.fHmpNumPhots), fHmpTrkFlags(c.fHmpTrkFlags),
 170   fN1(c.fN1),
 171   fN2(c.fN2),
 172   fPionEff(c.fPionEff),
 173   fKaonEff(c.fKaonEff),
 174   fProtEff(c.fProtEff),
 175   fPionTot(c.fPionTot),
 176   fKaonTot(c.fKaonTot),
 177   fProtTot(c.fProtTot),
 178   fPionNot(c.fPionNot),
 179   fKaonNot(c.fKaonNot),
 180   fProtNot(c.fProtNot),
 181   fPionCon(c.fPionCon),
 182   fKaonCon(c.fKaonCon),
 183   fProtCon(c.fProtCon),
 184   fThetavsPiFromK(c.fThetavsPiFromK),
 185   fThetapivsPesd(c.fThetapivsPesd),
 186   fThetaKvsPesd(c.fThetaKvsPesd),
 187   fThetaPvsPesd(c.fThetaPvsPesd)
 188 {
 189   //
 190   // Copy Constructor
 191   //
 192 }
 193
 194 //___________________________________________________________________________
 195 AliHMPIDAnalysisTask::~AliHMPIDAnalysisTask() {
 196   //
 197   //destructor
 198   //
 199   Info("~AliHMPIDAnalysisTask","Calling Destructor");
 200   if (fHmpHistList) {fHmpHistList->Clear(); delete fHmpHistList;}
 201 }
 202
 203 //___________________________________________________________________________
 204 void AliHMPIDAnalysisTask::ConnectInputData(Option_t *)
 205 {
 206   // Connect ESD here
 207
 208   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 209   if (!esdH) {
 210       AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
 211     } else
 212       fESD = esdH->GetEvent();
 213
 214   // Connect MC
 215
 216   AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 217   if (!mcH) {
 218       AliDebug(2,Form("ERROR: Could not get MCEventHandler"));
 219     } else
 220       fMC = mcH->MCEvent();
 221 }
 222
 223 //_________________________________________________
 224 void AliHMPIDAnalysisTask::Exec(Option_t *)
 225 {
 226   AliStack* pStack = fMC->Stack();
 227   Int_t label;
 228   Double_t priors[5]={1.,1.,1.,1.,1.}; //{0.01,0.01,0.83,0.10,0.5};
 229   Double_t probs[5];
 230   AliPID *pPid = new AliPID();
 231   pPid->SetPriors(priors);
 232   Double_t n = 1.293;
 233   Double_t dGeVMass[] = {0.000511,0.105658,0.13957018,0.493677,0.938272};
 234
 235   AliESDtrack *track=0;
 236   TParticle *pPart=0;
 237   //
 238   // Main loop function, executed on Event basis
 239   //
 240   for (Int_t iTrack = 0; iTrack<fESD->GetNumberOfTracks(); iTrack++) {
 241
 242     track = fESD->GetTrack(iTrack);
 243     if(!track) continue;
 244     Double_t dPtr = track->P();
 245     Double_t rin[3], rout[3];
 246     track->GetInnerXYZ(rin);
 247     track->GetOuterXYZ(rout);
 248     Double_t ktol = 0.001;
 249
 250     if(Equal(track->GetHMPIDsignal(),-20.,ktol))      fHmpTrkFlags->Fill(0);
 251     else if(Equal(track->GetHMPIDsignal(),-9.,ktol))  fHmpTrkFlags->Fill(1);
 252     else if(Equal(track->GetHMPIDsignal(),-5.,ktol))   fHmpTrkFlags->Fill(2);
 253     else if(Equal(track->GetHMPIDsignal(),-11.,ktol))  fHmpTrkFlags->Fill(3);
 254     else if(Equal(track->GetHMPIDsignal(),-22.,ktol)) fHmpTrkFlags->Fill(4);
 255     else fHmpTrkFlags->Fill(4);
 256
 257     if(Equal(track->GetHMPIDsignal(),-20.,ktol)) continue;
 258     if(track->GetHMPIDcluIdx() < 0) continue;
 259
 260     //Int_t ch = track->GetHMPIDcluIdx()/1000000;
 261     Int_t q, nph;
 262     Float_t x, y;
 263     Float_t xpc, ypc, th, ph;
 264     track->GetHMPIDmip(x,y,q,nph);
 265     track->GetHMPIDtrk(xpc,ypc,th,ph);
 266
 267     if(Equal(x,0.,ktol) && Equal(y,0.,ktol) && Equal(xpc,0.,ktol) && Equal(ypc,0.,ktol)) continue;
 268
 269     Double_t dist = TMath::Sqrt( (xpc-x)*(xpc-x) + (ypc - y)*(ypc - y));
 270     fHmpMipTrkDist->Fill(dist);
 271     fHmpMipTrkDistX->Fill(xpc - x);
 272     fHmpMipTrkDistY->Fill(ypc - y);
 273     Double_t pHmp[3] = {0}, pHmp3 = 0;
 274     if (track->GetOuterHmpPxPyPz(pHmp)) pHmp3 = TMath::Sqrt(pHmp[0]*pHmp[0]+pHmp[1]*pHmp[1]+pHmp[2]*pHmp[2]);
 275     if (dist <= 3.0) fHmpMipCharge3cm->Fill(q);
 276
 277     Float_t thetaTheorPion = 999., thetaTheorKaon = 999., thetaTheorProt = 999.,
 278             thetaHmpTheorPion = 999., thetaHmpTheorKaon = 999., thetaHmpTheorProt = 999.;
 279     if(dPtr != 0){
 280       thetaTheorPion = TMath::ACos(TMath::Sqrt(dPtr*dPtr + dGeVMass[2]*dGeVMass[2])/(n*dPtr));
 281       thetaTheorKaon = TMath::ACos(TMath::Sqrt(dPtr*dPtr + dGeVMass[3]*dGeVMass[3])/(n*dPtr));
 282       thetaTheorProt = TMath::ACos(TMath::Sqrt(dPtr*dPtr + dGeVMass[4]*dGeVMass[4])/(n*dPtr));
 283     }
 284     if(pHmp3 != 0){
 285       thetaHmpTheorPion = TMath::ACos(TMath::Sqrt(pHmp3*pHmp3 + dGeVMass[2]*dGeVMass[2])/(n*pHmp3));
 286       thetaHmpTheorKaon = TMath::ACos(TMath::Sqrt(pHmp3*pHmp3 + dGeVMass[3]*dGeVMass[3])/(n*pHmp3));
 287       thetaHmpTheorProt = TMath::ACos(TMath::Sqrt(pHmp3*pHmp3 + dGeVMass[4]*dGeVMass[4])/(n*pHmp3));
 288     }
 289
 290     track->GetHMPIDpid(probs);
 291     pPid->SetProbabilities(probs);
 292     if ((label = track->GetLabel()) < 0) continue;
 293     pPart = pStack->Particle(label);
 294
 295
 296     if(track->GetHMPIDsignal() > 0 )
 297     {
 298       fHmpPesdPhmp->Fill(track->P(),pHmp3);
 299       if (dist<=1.0) fHmpMipCharge1cm->Fill(q);
 300       fHmpNumPhots->Fill(nph);
 301       fHmpCkovPesd->Fill(track->P(),track->GetHMPIDsignal());
 302       fHmpCkovPhmp->Fill(pHmp3,track->GetHMPIDsignal());
 303
 304       if (!pStack->IsPhysicalPrimary(label)) continue;
 305       Int_t pdgCode = TMath::Abs(pPart->GetPdgCode());
 306       if (pdgCode==211){
 307         fThetapivsPesd->Fill(track->P(),track->GetHMPIDsignal());
 308         Int_t mot=pPart->GetFirstMother();
 309         if (mot > -1){
 310           TParticle *pMot=pStack->Particle(mot);
 311           TString str=pMot->GetName();
 312           if (str.Contains("K")) fThetavsPiFromK->Fill(pHmp3,track->GetHMPIDsignal());
 313         }
 314       }
 315       if (pdgCode==321) fThetaKvsPesd->Fill(track->P(),track->GetHMPIDsignal());
 316       if (pdgCode==2212) fThetaPvsPesd->Fill(track->P(),track->GetHMPIDsignal());
 317
 318       if (track->Pt()<1. || track->Pt()>5.) continue;
 319       Int_t ptBin=(Int_t) (2*(track->Pt()-1));
 320       if (pdgCode!=2212) fProtCon->Fill(ptBin);
 321       if (pdgCode==2212){
 322         fProtTot->Fill(ptBin);
 323         fProtEff->Fill(ptBin,pPid->GetProbability(AliPID::kProton));
 324         fPionNot->Fill(ptBin,pPid->GetProbability(AliPID::kPion));
 325         fKaonNot->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
 326       }
 327       if (pdgCode!=211) fPionCon->Fill(ptBin);
 328       if (pdgCode!=321) fKaonCon->Fill(ptBin);
 329       if (pdgCode==211){
 330         if (ptBin < 6){
 331           Float_t weight=pPid->GetProbability(AliPID::kElectron)+
 332                          pPid->GetProbability(AliPID::kMuon)+
 333                          pPid->GetProbability(AliPID::kPion);
 334           fPionTot->Fill(ptBin);
 335           fPionEff->Fill(ptBin,weight);
 336           fKaonNot->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
 337         }
 338         fProtNot->Fill(ptBin,pPid->GetProbability(AliPID::kProton));
 339       }
 340       if (pdgCode==321){
 341         if (ptBin < 6){
 342           fKaonTot->Fill(ptBin);
 343           fKaonEff->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
 344           fPionNot->Fill(ptBin,(pPid->GetProbability(AliPID::kPion)));
 345         }
 346         fProtNot->Fill(ptBin,(pPid->GetProbability(AliPID::kProton)));
 347       }
 348     }//there is signal
 349   }//track loop
 350
 351   delete pPid;
 352
 353   /* PostData(0) is taken care of by AliAnalysisTaskSE */
 354   PostData(0,fHmpHistList);
 355 }
 356
 357 //___________________________________________________________________________
 358 void AliHMPIDAnalysisTask::Terminate(Option_t*)
 359 {
 360   // The Terminate() function is the last function to be called during
 361   // a query. It always runs on the client, it can be used to present
 362   // the results graphically or save the results to file.
 363
 364   Info("Terminate","");
 365   fHmpHistList = dynamic_cast<TList*> (GetOutputData(0));
 366
 367   fPionEff = dynamic_cast<TH1F*> (fHmpHistList->FindObject("PionEff"));
 368   fPionTot = dynamic_cast<TH1I*> (fHmpHistList->FindObject("PionTot"));
 369   fPionNot = dynamic_cast<TH1F*> (fHmpHistList->FindObject("PionNot"));
 370   fPionCon = dynamic_cast<TH1I*> (fHmpHistList->FindObject("PionCon"));
 371   fKaonEff = dynamic_cast<TH1F*> (fHmpHistList->FindObject("KaonEff"));
 372   fKaonTot = dynamic_cast<TH1I*> (fHmpHistList->FindObject("KaonTot"));
 373   fKaonNot = dynamic_cast<TH1F*> (fHmpHistList->FindObject("KaonNot"));
 374   fKaonCon = dynamic_cast<TH1I*> (fHmpHistList->FindObject("KaonCon"));
 375   fProtEff = dynamic_cast<TH1F*> (fHmpHistList->FindObject("ProtEff"));
 376   fProtTot = dynamic_cast<TH1I*> (fHmpHistList->FindObject("ProtTot"));
 377   fProtNot = dynamic_cast<TH1F*> (fHmpHistList->FindObject("ProtNot"));
 378   fProtCon = dynamic_cast<TH1I*> (fHmpHistList->FindObject("ProtCon"));
 379
 380   Float_t *pionEff=fPionEff->GetArray();
 381   Int_t   *pionTot=fPionTot->GetArray();
 382   Float_t *pionNot=fPionNot->GetArray();
 383   Int_t   *pionCon=fPionCon->GetArray();
 384   Float_t *kaonEff=fKaonEff->GetArray();
 385   Int_t   *kaonTot=fKaonTot->GetArray();
 386   Float_t *kaonNot=fKaonNot->GetArray();
 387   Int_t   *kaonCon=fKaonCon->GetArray();
 388   Float_t *protEff=fProtEff->GetArray();
 389   Int_t   *protTot=fProtTot->GetArray();
 390   Float_t *protNot=fProtNot->GetArray();
 391   Int_t   *protCon=fProtCon->GetArray();
 392
 393   TGraphErrors *effPi = new TGraphErrors(fN1);
 394   TGraphErrors *effKa = new TGraphErrors(fN1);
 395   TGraphErrors *effPr = new TGraphErrors(fN2);
 396   TGraphErrors *conPi = new TGraphErrors(fN1);
 397   TGraphErrors *conKa = new TGraphErrors(fN1);
 398   TGraphErrors *conPr = new TGraphErrors(fN2);
 399
 400   Float_t pt[8]={1.25,1.75,2.25,2.75,3.25,3.75,4.25,4.75};
 401   Float_t eff=0, effErr=0, con=0, conErr=0;
 402   for (Int_t i=0; i< fN2; i++){
 403     eff = protEff[i+1]/TMath::Max(protTot[i+1],1);
 404     effErr = TMath::Sqrt(eff*(1.-eff)/TMath::Max(protTot[i+1],1));
 405     con = protNot[i+1]/TMath::Max(protCon[i+1],1);
 406     conErr = TMath::Sqrt(con*(1.-con)/protCon[i+1]);
 407     effPr->SetPoint(i,pt[i],eff);
 408     effPr->SetPointError(i,0,effErr);
 409     conPr->SetPoint(i,pt[i],con);
 410     conPr->SetPointError(i,0,conErr);
 411
 412     if (i>=fN1) continue;
 413     eff = pionEff[i+1]/pionTot[i+1];
 414     effErr = TMath::Sqrt(eff*(1.-eff)/pionTot[i+1]);
 415     con=pionNot[i+1]/pionCon[i+1];
 416     conErr = TMath::Sqrt(con*(1.-con)/pionCon[i+1]);
 417     effPi->SetPoint(i,pt[i],(Float_t)pionEff[i+1]/(Float_t)pionTot[i+1]);
 418     effPi->SetPointError(i,0,effErr);
 419     conPi->SetPoint(i,pt[i],(Float_t)pionNot[i+1]/(Float_t)pionCon[i+1]);
 420     conPi->SetPointError(i,0,conErr);
 421
 422     eff = kaonEff[i+1]/TMath::Max(kaonTot[i+1],1);
 423     effErr = TMath::Sqrt(eff*(1.-eff)/kaonTot[i+1]);
 424     con = kaonNot[i+1]/TMath::Max(kaonCon[i+1],1);
 425     conErr = TMath::Sqrt(con*(1.-con)/kaonCon[i+1]);
 426     effKa->SetPoint(i,pt[i],(Float_t)kaonEff[i+1]/TMath::Max(kaonTot[i+1],1));
 427     effKa->SetPointError(i,0,effErr);
 428     conKa->SetPoint(i,pt[i],(Float_t)kaonNot[i+1]/TMath::Max(kaonCon[i+1],1));
 429     conKa->SetPointError(i,0,conErr);
 430   }
 431
 432   TCanvas *pCan=new TCanvas("Hmp","Efficiency and contamination",500,900);
 433   pCan->Divide(1,3);
 434
 435   pCan->cd(1);
 436   effPi->SetTitle("Pions");
 437   effPi->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 438   effPi->GetYaxis()->SetRangeUser(0.,1.);
 439   effPi->SetMarkerStyle(20);
 440   effPi->Draw("ALP");
 441   conPi->SetMarkerStyle(21);
 442   conPi->Draw("sameLP");
 443
 444   pCan->cd(2);
 445   effKa->SetTitle("Kaons");
 446   effKa->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 447   effKa->GetYaxis()->SetRangeUser(0.,1.);
 448   effKa->SetMarkerStyle(20);
 449   effKa->Draw("ALP");
 450   conKa->SetMarkerStyle(21);
 451   conKa->Draw("sameLP");
 452
 453   pCan->cd(3);
 454   effPr->SetTitle("Protons");
 455   effPr->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 456   effPr->GetYaxis()->SetRangeUser(0.,1.);
 457   effPr->SetMarkerStyle(20);
 458   effPr->Draw("ALP");
 459   conPr->SetMarkerStyle(21);
 460   conPr->Draw("sameLP");
 461
 462   TFile *outFile = new TFile("HmpidGraphs.root","recreate");
 463   pCan->Write();
 464   outFile->Close();
 465
 466   AliAnalysisTaskSE::Terminate();
 467
 468 }
 469
 470 //___________________________________________________________________________
 471 void AliHMPIDAnalysisTask::CreateOutputObjects() {
 472   //
 473   //HERE ONE CAN CREATE OUTPUT OBJECTS
 474   //TO BE SET BEFORE THE EXECUTION OF THE TASK
 475   //
 476
 477   //slot #1
 478    OpenFile(0);
 479    fHmpHistList = new TList();
 480
 481    fHmpPesdPhmp = new TH2F("fHmpPesdPhmp","HMPID: ESD p - running p;HMP p (GeV/c);ESD p (Gev/c)",100,0,10,100,0,10);
 482    fHmpHistList->Add(fHmpPesdPhmp);
 483
 484    fHmpCkovPesd = new TH2F("fHmpCkovPesd","HMPID: ThetaCherenkov vs P;p_esd (GeV/c);#Theta_C;Entries",100,0,10,110,0,1.1);
 485    fHmpHistList->Add(fHmpCkovPesd);
 486
 487    fHmpCkovPhmp = new TH2F("fHmpCkovPhmp","HMPID: ThetaCherenkov vs P;p_hmp (GeV/c);#Theta_C;Entries",100,0,10,110,0,1.1);
 488    fHmpHistList->Add(fHmpCkovPhmp);
 489
 490    fHmpMipTrkDist = new TH1F("fHmpMipTrkDist","HMPID MIP-Track distance;distance (cm);Entries",800,-20,20);
 491    fHmpHistList->Add(fHmpMipTrkDist);
 492
 493    fHmpMipTrkDistX = new TH1F("fHmpMipTrkDistX","HMPID MIP-Track distance in local X;distance (cm);Entries",800,-20,20);
 494    fHmpHistList->Add(fHmpMipTrkDistX);
 495
 496    fHmpMipTrkDistY = new TH1F("fHmpMipTrkDistY","HMPID MIP-Track distance in local Y;distance (cm);Entries",800,-20,20);
 497    fHmpHistList->Add(fHmpMipTrkDistY);
 498
 499    fHmpMipCharge3cm = new TH1F("fHmpMipCharge3cm","HMPID MIP Charge;MIP Charge (ADC);Entries",5001,-0.5,5000.5);
 500    fHmpHistList->Add(fHmpMipCharge3cm);
 501
 502    fHmpMipCharge1cm = new TH1F("fHmpMipCharge1cm","HMPID MIP Charge;MIP Charge (ADC);Entries",5001,-0.5,5000.5);
 503    fHmpHistList->Add(fHmpMipCharge1cm);
 504
 505    fHmpNumPhots = new TH1F("fHmpNumPhots","HMPID Number of photon clusters on ring;#photon clus.;Entries",51,-0.5,50.5);
 506    fHmpHistList->Add(fHmpNumPhots);
 507
 508    fHmpTrkFlags = new TH1F("fHmpTrkFlags","HMPID track flags",6,0,6);
 509    TString summary[6] =  {"NotPerformed","MipDistCut", "MipQdcCut", "NoPhotAccept", "kNoRad", "other"};
 510    for(Int_t ibin = 0; ibin < 6; ibin++) fHmpTrkFlags->GetXaxis()->SetBinLabel(ibin+1,Form("%i  %s",ibin+1,summary[ibin].Data()));
 511    fHmpHistList->Add(fHmpTrkFlags);
 512
 513    fPionEff = new TH1F("PionEff","Identified pions",fN1,0,fN1);
 514    fKaonEff = new TH1F("KaonEff","Identified kaons",fN1,0,fN1);
 515    fProtEff = new TH1F("ProtEff","Identified protons",fN2,0,fN2);
 516    fPionTot = new TH1I("PionTot","Total MC pions",fN1,0,fN1);
 517    fKaonTot = new TH1I("KaonTot","Total MC kaons",fN1,0,fN1);
 518    fProtTot = new TH1I("ProtTot","Total MC protons",fN2,0,fN2);
 519    fPionNot = new TH1F("PionNot","Misidentified pions",fN1,0,fN1);
 520    fKaonNot = new TH1F("KaonNot","Misidentified kaons",fN1,0,fN1);
 521    fProtNot = new TH1F("ProtNot","Misidentified protons",fN2,0,fN2);
 522    fPionCon = new TH1I("PionCon","Total not MC pions",fN1,0,fN1);
 523    fKaonCon = new TH1I("KaonCon","Total not MC kaons",fN1,0,fN1);
 524    fProtCon = new TH1I("ProtCon","Total not MC protons",fN2,0,fN2);
 525
 526    fHmpHistList->Add(fPionEff); fHmpHistList->Add(fKaonEff); fHmpHistList->Add(fProtEff);
 527    fHmpHistList->Add(fPionTot); fHmpHistList->Add(fKaonTot); fHmpHistList->Add(fProtTot);
 528    fHmpHistList->Add(fPionNot); fHmpHistList->Add(fKaonNot); fHmpHistList->Add(fProtNot);
 529    fHmpHistList->Add(fPionCon); fHmpHistList->Add(fKaonCon); fHmpHistList->Add(fProtCon);
 530
 531    fThetavsPiFromK= new TH2F("ThetavsPiFromK","Theta vs p of pions from K;p_esd (GeV/c);#Theta_C",100,0,10,110,0,1.1);
 532    fHmpHistList->Add(fThetavsPiFromK);
 533
 534    fThetapivsPesd = new TH2F("ThetapivsPesd","Theta of pions vs p of esd;p_esd (GeV/c);#Theta_C",100,0,10,110,0,1.1);
 535    fHmpHistList->Add(fThetapivsPesd);
 536
 537    fThetaKvsPesd  = new TH2F("ThetaKvsPesd","Theta of kaons vs p of esd;p_esd (GeV/c);#Theta_C",100,0,10,110,0,1.1);
 538    fHmpHistList->Add(fThetaKvsPesd);
 539
 540    fThetaPvsPesd  = new TH2F("ThetaPvsPesd","Theta of protons vs p of esd;p_esd (GeV/c);#Theta_C",100,0,10,110,0,1.1);
 541    fHmpHistList->Add(fThetaPvsPesd);
 542
 543 }
 544 //____________________________________________________________________________________________________________________________________
 545 Bool_t AliHMPIDAnalysisTask::Equal(Double_t x, Double_t y, Double_t tolerance)
 546 {
 547  return abs(x - y) <= tolerance ;
 548 }
 549
 550 #endif