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
  18 // A set of histograms is created.
  19 //==============================================================================
  20 //
  21 // By means of AliHMPIDAnalysisTask.C macro it is possible to use this class
  22 // to perform the analysis on local data, on data on alien using local machine
  23 // and on CAF.
  24
  25 #ifndef AliHMPIDAnalysisTASK_CXX
  26 #define AliHMPIDAnalysisTASK_CXX
  27
  28
  29 #include "TH1.h"
  30 #include "TH2.h"
  31 #include "TFile.h"
  32 #include "TCanvas.h"
  33 #include "TGraphErrors.h"
  34 #include "AliAnalysisManager.h"
  35 #include "AliESDInputHandler.h"
  36 #include "AliMCEventHandler.h"
  37 #include "AliMCEvent.h"
  38 #include "AliESDtrack.h"
  39 #include "AliPID.h"
  40 #include "AliLog.h"
  41 #include "AliHMPIDAnalysisTask.h"
  42
  43 ClassImp(AliHMPIDAnalysisTask)
  44
  45 //__________________________________________________________________________
  46 AliHMPIDAnalysisTask::AliHMPIDAnalysisTask() :
  47   fESD(0x0),fMC(0x0),fUseMC(kTRUE),
  48   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   fTree(0x0)
  74 {
  75   //
  76   //Default ctor
  77   //
  78   for (Int_t i=0; i<34; i++) fVar[i]=0;
  79 }
  80
  81 //___________________________________________________________________________
  82 AliHMPIDAnalysisTask::AliHMPIDAnalysisTask(const Char_t* name) :
  83   AliAnalysisTaskSE(name),
  84   fESD(0x0), fMC(0x0), fUseMC(kTRUE),
  85   fHmpHistList(0x0), fNevts(0),
  86   fTrigNevts(0),
  87   fTrigger(0),
  88   fHmpPesdPhmp(0x0), fHmpCkovPesd(0x0), fHmpCkovPhmp(0x0),
  89   fHmpMipTrkDist(0x0), fHmpMipTrkDistX(0x0), fHmpMipTrkDistY(0x0),
  90   fHmpMipCharge3cm(0x0), fHmpMipCharge1cm(0x0),
  91   fHmpNumPhots(0x0), fHmpTrkFlags(0x0),
  92   fN1(6),
  93   fN2(8),
  94   fPionEff(0x0),
  95   fKaonEff(0x0),
  96   fProtEff(0x0),
  97   fPionTot(0x0),
  98   fKaonTot(0x0),
  99   fProtTot(0x0),
 100   fPionNot(0x0),
 101   fKaonNot(0x0),
 102   fProtNot(0x0),
 103   fPionCon(0x0),
 104   fKaonCon(0x0),
 105   fProtCon(0x0),
 106   fThetavsPiFromK(0x0),
 107   fThetapivsPesd(0x0),
 108   fThetaKvsPesd(0x0),
 109   fThetaPvsPesd(0x0),
 110   fTree(0x0)
 111 {
 112   //
 113   // Constructor. Initialization of Inputs and Outputs
 114   //
 115   for (Int_t i=0; i<34; i++) fVar[i]=0;
 116
 117   DefineOutput(1,TList::Class());
 118   DefineOutput(2,TTree::Class());
 119 }
 120
 121 //___________________________________________________________________________
 122 AliHMPIDAnalysisTask& AliHMPIDAnalysisTask::operator=(const AliHMPIDAnalysisTask& c)
 123 {
 124   //
 125   // Assignment operator
 126   //
 127   if (this!=&c) {
 128     AliAnalysisTaskSE::operator=(c);
 129     fESD             = c.fESD;
 130     fMC              = c.fMC;
 131     fUseMC           = c.fUseMC;
 132     fHmpHistList     = c.fHmpHistList;
 133     fNevts           = c.fNevts;
 134     fTrigNevts       = c.fTrigNevts;
 135     fTrigger         = c.fTrigger;
 136     fHmpPesdPhmp     = c.fHmpPesdPhmp;
 137     fHmpCkovPesd     = c.fHmpCkovPesd;
 138     fHmpCkovPhmp     = c.fHmpCkovPhmp;
 139     fHmpMipTrkDist   = c.fHmpMipTrkDist;
 140     fHmpMipTrkDistX  = c.fHmpMipTrkDistX;
 141     fHmpMipTrkDistY  = c.fHmpMipTrkDistY;
 142     fHmpMipCharge3cm = c.fHmpMipCharge3cm;
 143     fHmpMipCharge1cm = c.fHmpMipCharge1cm;
 144     fHmpNumPhots     = c.fHmpNumPhots;
 145     fHmpTrkFlags     = c.fHmpTrkFlags;
 146     fN1              = c.fN1;
 147     fN2              = c.fN2;
 148     fPionEff         = c.fPionEff;
 149     fKaonEff         = c.fKaonEff;
 150     fProtEff         = c.fProtEff;
 151     fPionTot         = c.fPionTot;
 152     fKaonTot         = c.fKaonTot;
 153     fProtTot         = c.fProtTot;
 154     fPionNot         = c.fPionNot;
 155     fKaonNot         = c.fKaonNot;
 156     fProtNot         = c.fProtNot;
 157     fPionCon         = c.fPionCon;
 158     fKaonCon         = c.fKaonCon;
 159     fProtCon         = c.fProtCon;
 160     fThetavsPiFromK  = c.fThetavsPiFromK;
 161     fThetapivsPesd   = c.fThetapivsPesd;
 162     fThetaKvsPesd    = c.fThetaKvsPesd;
 163     fThetaPvsPesd    = c.fThetaPvsPesd;
 164     fTree            = c.fTree;
 165     for (Int_t i=0; i<34; i++) fVar[i]=c.fVar[i];
 166   }
 167   return *this;
 168 }
 169
 170 //___________________________________________________________________________
 171 AliHMPIDAnalysisTask::AliHMPIDAnalysisTask(const AliHMPIDAnalysisTask& c) :
 172   AliAnalysisTaskSE(c),
 173   fESD(c.fESD),fMC(c.fMC),fUseMC(c.fUseMC),
 174   fHmpHistList(c.fHmpHistList), fNevts(c.fNevts),
 175   fTrigNevts(c.fTrigNevts),
 176   fTrigger(c.fTrigger),
 177   fHmpPesdPhmp(c.fHmpPesdPhmp),fHmpCkovPesd(c.fHmpCkovPesd),fHmpCkovPhmp(c.fHmpCkovPhmp),
 178   fHmpMipTrkDist(c.fHmpMipTrkDist),fHmpMipTrkDistX(c.fHmpMipTrkDistX),fHmpMipTrkDistY(c.fHmpMipTrkDistY),
 179   fHmpMipCharge3cm(c.fHmpMipCharge3cm), fHmpMipCharge1cm(c.fHmpMipCharge1cm),
 180   fHmpNumPhots(c.fHmpNumPhots), fHmpTrkFlags(c.fHmpTrkFlags),
 181   fN1(c.fN1),
 182   fN2(c.fN2),
 183   fPionEff(c.fPionEff),
 184   fKaonEff(c.fKaonEff),
 185   fProtEff(c.fProtEff),
 186   fPionTot(c.fPionTot),
 187   fKaonTot(c.fKaonTot),
 188   fProtTot(c.fProtTot),
 189   fPionNot(c.fPionNot),
 190   fKaonNot(c.fKaonNot),
 191   fProtNot(c.fProtNot),
 192   fPionCon(c.fPionCon),
 193   fKaonCon(c.fKaonCon),
 194   fProtCon(c.fProtCon),
 195   fThetavsPiFromK(c.fThetavsPiFromK),
 196   fThetapivsPesd(c.fThetapivsPesd),
 197   fThetaKvsPesd(c.fThetaKvsPesd),
 198   fThetaPvsPesd(c.fThetaPvsPesd),
 199   fTree(c.fTree)
 200 {
 201   //
 202   // Copy Constructor
 203   //
 204   for (Int_t i=0; i<34; i++) fVar[i]=c.fVar[i];
 205 }
 206
 207 //___________________________________________________________________________
 208 AliHMPIDAnalysisTask::~AliHMPIDAnalysisTask() {
 209   //
 210   //destructor
 211   //
 212   Info("~AliHMPIDAnalysisTask","Calling Destructor");
 213   if (fHmpHistList /*&& !AliAnalysisManager::GetAnalysisManager()->IsProofMode()*/) delete fHmpHistList;
 214 }
 215
 216 //___________________________________________________________________________
 217 void AliHMPIDAnalysisTask::ConnectInputData(Option_t *)
 218 {
 219   // Connect ESD here
 220
 221   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 222   if (!esdH) {
 223     AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
 224   } else
 225     fESD = esdH->GetEvent();
 226
 227   if (fUseMC){
 228     // Connect MC
 229     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 230     if (!mcH) {
 231       AliDebug(2,Form("ERROR: Could not get MCEventHandler"));
 232       fUseMC = kFALSE;
 233     } else
 234       fMC = mcH->MCEvent();
 235       if (!fMC) AliDebug(2,Form("ERROR: Could not get MCEvent"));
 236   }
 237 }
 238
 239 //___________________________________________________________________________
 240 void AliHMPIDAnalysisTask::UserExec(Option_t *)
 241 {
 242   Double_t priors[5]={1.,1.,1.,1.,1.}; //{0.01,0.01,0.83,0.10,0.5};
 243   Double_t probs[5];
 244   AliPID *pPid = new AliPID();
 245   pPid->SetPriors(priors);
 246   Double_t n = 1.293;
 247   Double_t dGeVMass[] = {0.000511,0.105658,0.13957018,0.493677,0.938272};
 248   AliESDtrack *track=0;
 249   TParticle *pPart=0;
 250   AliStack* pStack = 0;
 251   Int_t label;
 252   if (fUseMC){
 253     pStack = fMC->Stack();
 254   }
 255
 256   //
 257   // Main loop function, executed on Event basis
 258   //
 259   for (Int_t iTrack = 0; iTrack<fESD->GetNumberOfTracks(); iTrack++) {
 260
 261     track = fESD->GetTrack(iTrack);
 262     if(!track) continue;
 263     Double_t dPtr = track->P();
 264     Double_t rin[3], rout[3];
 265     track->GetInnerXYZ(rin);
 266     track->GetOuterXYZ(rout);
 267     Double_t ktol = 0.001;
 268
 269     if(Equal(track->GetHMPIDsignal(),-20.,ktol))      fHmpTrkFlags->Fill(0);
 270     else if(Equal(track->GetHMPIDsignal(),-9.,ktol))  fHmpTrkFlags->Fill(1);
 271     else if(Equal(track->GetHMPIDsignal(),-5.,ktol))  fHmpTrkFlags->Fill(2);
 272     else if(Equal(track->GetHMPIDsignal(),-11.,ktol)) fHmpTrkFlags->Fill(3);
 273     else if(Equal(track->GetHMPIDsignal(),-22.,ktol)) fHmpTrkFlags->Fill(4);
 274     else fHmpTrkFlags->Fill(4);
 275
 276     if(Equal(track->GetHMPIDsignal(),-20.,ktol)) continue;
 277     if(track->GetHMPIDcluIdx() < 0) continue;
 278
 279     Int_t q, nph;
 280     Float_t x, y;
 281     Float_t xpc, ypc, th, ph;
 282     track->GetHMPIDmip(x,y,q,nph);
 283     track->GetHMPIDtrk(xpc,ypc,th,ph);
 284
 285     if(Equal(x,0.,ktol) && Equal(y,0.,ktol) && Equal(xpc,0.,ktol) && Equal(ypc,0.,ktol)) continue;
 286
 287     Double_t dist = TMath::Sqrt( (xpc-x)*(xpc-x) + (ypc - y)*(ypc - y));
 288     fHmpMipTrkDist->Fill(dist);
 289     fHmpMipTrkDistX->Fill(xpc - x);
 290     fHmpMipTrkDistY->Fill(ypc - y);
 291     Double_t pHmp[3] = {0}, pHmp3 = 0;
 292     if (track->GetOuterHmpPxPyPz(pHmp)) pHmp3 = TMath::Sqrt(pHmp[0]*pHmp[0]+pHmp[1]*pHmp[1]+pHmp[2]*pHmp[2]);
 293     if (dist <= 3.0) fHmpMipCharge3cm->Fill(q);
 294
 295     Float_t thetaTheorPion = 999., thetaTheorKaon = 999., thetaTheorProt = 999.,
 296             thetaHmpTheorPion = 999., thetaHmpTheorKaon = 999., thetaHmpTheorProt = 999.;
 297     if(dPtr != 0){
 298       thetaTheorPion = TMath::ACos(TMath::Sqrt(dPtr*dPtr + dGeVMass[2]*dGeVMass[2])/(n*dPtr));
 299       thetaTheorKaon = TMath::ACos(TMath::Sqrt(dPtr*dPtr + dGeVMass[3]*dGeVMass[3])/(n*dPtr));
 300       thetaTheorProt = TMath::ACos(TMath::Sqrt(dPtr*dPtr + dGeVMass[4]*dGeVMass[4])/(n*dPtr));
 301     }
 302     if(pHmp3 != 0){
 303       thetaHmpTheorPion = TMath::ACos(TMath::Sqrt(pHmp3*pHmp3 + dGeVMass[2]*dGeVMass[2])/(n*pHmp3));
 304       thetaHmpTheorKaon = TMath::ACos(TMath::Sqrt(pHmp3*pHmp3 + dGeVMass[3]*dGeVMass[3])/(n*pHmp3));
 305       thetaHmpTheorProt = TMath::ACos(TMath::Sqrt(pHmp3*pHmp3 + dGeVMass[4]*dGeVMass[4])/(n*pHmp3));
 306     }
 307
 308     Float_t b[2];
 309     Float_t bCov[3];
 310     track->GetImpactParameters(b,bCov);
 311
 312     track->GetHMPIDpid(probs);
 313     pPid->SetProbabilities(probs);
 314     if (fUseMC){
 315       if ((label = track->GetLabel()) < 0) continue;
 316       pPart = pStack->Particle(label);
 317     }
 318
 319     if(track->GetHMPIDsignal() > 0 ){
 320       if (pHmp3) fHmpPesdPhmp->Fill(track->P(),pHmp3);
 321       if (dist<=1.0) fHmpMipCharge1cm->Fill(q);
 322       fHmpNumPhots->Fill(nph);
 323       fHmpCkovPesd->Fill(track->P(),track->GetHMPIDsignal());
 324       if (pHmp3) fHmpCkovPhmp->Fill(pHmp3,track->GetHMPIDsignal());
 325
 326       if (fUseMC && dist<0.5 && TMath::Abs(th)<0.13){
 327         if (!pStack->IsPhysicalPrimary(label)) continue;
 328         Int_t pdgCode = TMath::Abs(pPart->GetPdgCode());
 329         if (pdgCode==211){
 330           fThetapivsPesd->Fill(track->P(),track->GetHMPIDsignal());
 331           Int_t mot=pPart->GetFirstMother();
 332           if (mot > -1){
 333             TParticle *pMot=pStack->Particle(mot);
 334             TString str=pMot->GetName();
 335             if (str.Contains("K")) fThetavsPiFromK->Fill(pHmp3,track->GetHMPIDsignal());
 336           }
 337         }
 338         if (pdgCode==321) fThetaKvsPesd->Fill(track->P(),track->GetHMPIDsignal());
 339         if (pdgCode==2212) fThetaPvsPesd->Fill(track->P(),track->GetHMPIDsignal());
 340
 341         if (track->Pt()<1. || track->Pt()>5.) continue;
 342         Int_t ptBin=(Int_t) (2*(track->Pt()-1));
 343         if (pdgCode!=2212) fProtCon->Fill(ptBin);
 344         if (pdgCode==2212){
 345           fProtTot->Fill(ptBin);
 346           fProtEff->Fill(ptBin,pPid->GetProbability(AliPID::kProton));
 347           fPionNot->Fill(ptBin,pPid->GetProbability(AliPID::kPion));
 348           fKaonNot->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
 349         }
 350         if (pdgCode!=211) fPionCon->Fill(ptBin);
 351         if (pdgCode!=321) fKaonCon->Fill(ptBin);
 352         if (pdgCode==211){
 353           if (ptBin < 6){
 354             Float_t weight=pPid->GetProbability(AliPID::kElectron)+
 355                            pPid->GetProbability(AliPID::kMuon)+
 356                            pPid->GetProbability(AliPID::kPion);
 357             fPionTot->Fill(ptBin);
 358             fPionEff->Fill(ptBin,weight);
 359             fKaonNot->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
 360           }
 361           fProtNot->Fill(ptBin,pPid->GetProbability(AliPID::kProton));
 362         }
 363         if (pdgCode==321){
 364           if (ptBin < 6){
 365             fKaonTot->Fill(ptBin);
 366             fKaonEff->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
 367             fPionNot->Fill(ptBin,(pPid->GetProbability(AliPID::kPion)));
 368           }
 369           fProtNot->Fill(ptBin,(pPid->GetProbability(AliPID::kProton)));
 370         }
 371       }
 372     }//there is signal
 373     fVar[0] = track->GetHMPIDcluIdx()/1000000;
 374     fVar[1] = pHmp3;
 375     fVar[2] = dPtr;
 376     fVar[3] = xpc;
 377     fVar[4] = ypc;
 378     fVar[5] = x;
 379     fVar[6] = y;
 380     fVar[7] = thetaTheorPion;
 381     fVar[8] = thetaTheorKaon;
 382     fVar[9] = thetaTheorProt;
 383     fVar[10] = thetaHmpTheorPion;
 384     fVar[11] = thetaHmpTheorKaon;
 385     fVar[12] = thetaHmpTheorProt;
 386     fVar[13] = (Float_t)track->GetHMPIDsignal();
 387     fVar[14] = q;
 388     fVar[15] = th;
 389     fVar[16] = ph;
 390     fVar[17] = (Float_t)track->GetSign();
 391     fVar[18] = (Float_t)nph;
 392     fVar[19] = (Float_t)track->GetNcls(1);
 393     fVar[20] = (Float_t)probs[0];
 394     fVar[21] = (Float_t)probs[1];
 395     fVar[22] = (Float_t)probs[2];
 396     fVar[23] = (Float_t)probs[3];
 397     fVar[24] = (Float_t)probs[4];
 398     fVar[25] = (Float_t)track->GetTOFsignal();
 399     fVar[26] = (Float_t)track->GetKinkIndex(0);
 400     fVar[27] = (Float_t)track->Xv();
 401     fVar[28] = (Float_t)track->Yv();
 402     fVar[29] = (Float_t)track->Zv();
 403     fVar[30] = (Float_t)track->GetTPCchi2();
 404     fVar[31] = b[0];
 405     fVar[32] = b[1];
 406     fVar[33] = track->GetHMPIDcluIdx()%1000000/1000;
 407     fTree->Fill();
 408   }//track loop
 409   delete pPid;
 410
 411   /* PostData(0) is taken care of by AliAnalysisTaskSE */
 412   PostData(1,fHmpHistList);
 413   PostData(2,fTree);
 414 }
 415
 416 //___________________________________________________________________________
 417 void AliHMPIDAnalysisTask::Terminate(Option_t*)
 418 {
 419   // The Terminate() function is the last function to be called during
 420   // a query. It always runs on the client, it can be used to present
 421   // the results graphically or save the results to file.
 422
 423   Info("Terminate","");
 424
 425   if (!fUseMC) return;
 426
 427   fHmpHistList = dynamic_cast<TList*> (GetOutputData(1));
 428
 429   if (!fHmpHistList) {
 430     AliError("Histogram List is not available");
 431     return;
 432   }
 433
 434   fPionEff = dynamic_cast<TH1F*> (fHmpHistList->FindObject("PionEff"));
 435   fPionTot = dynamic_cast<TH1I*> (fHmpHistList->FindObject("PionTot"));
 436   fPionNot = dynamic_cast<TH1F*> (fHmpHistList->FindObject("PionNot"));
 437   fPionCon = dynamic_cast<TH1I*> (fHmpHistList->FindObject("PionCon"));
 438   fKaonEff = dynamic_cast<TH1F*> (fHmpHistList->FindObject("KaonEff"));
 439   fKaonTot = dynamic_cast<TH1I*> (fHmpHistList->FindObject("KaonTot"));
 440   fKaonNot = dynamic_cast<TH1F*> (fHmpHistList->FindObject("KaonNot"));
 441   fKaonCon = dynamic_cast<TH1I*> (fHmpHistList->FindObject("KaonCon"));
 442   fProtEff = dynamic_cast<TH1F*> (fHmpHistList->FindObject("ProtEff"));
 443   fProtTot = dynamic_cast<TH1I*> (fHmpHistList->FindObject("ProtTot"));
 444   fProtNot = dynamic_cast<TH1F*> (fHmpHistList->FindObject("ProtNot"));
 445   fProtCon = dynamic_cast<TH1I*> (fHmpHistList->FindObject("ProtCon"));
 446
 447   Float_t *pionEff=fPionEff->GetArray();
 448   Int_t   *pionTot=fPionTot->GetArray();
 449   Float_t *pionNot=fPionNot->GetArray();
 450   Int_t   *pionCon=fPionCon->GetArray();
 451   Float_t *kaonEff=fKaonEff->GetArray();
 452   Int_t   *kaonTot=fKaonTot->GetArray();
 453   Float_t *kaonNot=fKaonNot->GetArray();
 454   Int_t   *kaonCon=fKaonCon->GetArray();
 455   Float_t *protEff=fProtEff->GetArray();
 456   Int_t   *protTot=fProtTot->GetArray();
 457   Float_t *protNot=fProtNot->GetArray();
 458   Int_t   *protCon=fProtCon->GetArray();
 459
 460   TGraphErrors *effPi = new TGraphErrors(fN1);
 461   TGraphErrors *effKa = new TGraphErrors(fN1);
 462   TGraphErrors *effPr = new TGraphErrors(fN2);
 463   TGraphErrors *conPi = new TGraphErrors(fN1);
 464   TGraphErrors *conKa = new TGraphErrors(fN1);
 465   TGraphErrors *conPr = new TGraphErrors(fN2);
 466
 467   Float_t pt[8]={1.25,1.75,2.25,2.75,3.25,3.75,4.25,4.75};
 468   Float_t eff=0, effErr=0, con=0, conErr=0;
 469   for (Int_t i=0; i< fN2; i++){
 470     eff = protEff[i+1]/TMath::Max(protTot[i+1],1);
 471     effErr = TMath::Sqrt(eff*(1.-eff)/TMath::Max(protTot[i+1],1));
 472     con = protNot[i+1]/TMath::Max(protCon[i+1],1);
 473     conErr = TMath::Sqrt(con*(1.-con)/protCon[i+1]);
 474     effPr->SetPoint(i,pt[i],eff);
 475     effPr->SetPointError(i,0,effErr);
 476     conPr->SetPoint(i,pt[i],con);
 477     conPr->SetPointError(i,0,conErr);
 478
 479     if (i>=fN1) continue;
 480     eff = pionEff[i+1]/pionTot[i+1];
 481     effErr = TMath::Sqrt(eff*(1.-eff)/pionTot[i+1]);
 482     con=pionNot[i+1]/pionCon[i+1];
 483     conErr = TMath::Sqrt(con*(1.-con)/pionCon[i+1]);
 484     effPi->SetPoint(i,pt[i],(Float_t)pionEff[i+1]/(Float_t)pionTot[i+1]);
 485     effPi->SetPointError(i,0,effErr);
 486     conPi->SetPoint(i,pt[i],(Float_t)pionNot[i+1]/(Float_t)pionCon[i+1]);
 487     conPi->SetPointError(i,0,conErr);
 488
 489     eff = kaonEff[i+1]/TMath::Max(kaonTot[i+1],1);
 490     effErr = TMath::Sqrt(eff*(1.-eff)/kaonTot[i+1]);
 491     con = kaonNot[i+1]/TMath::Max(kaonCon[i+1],1);
 492     conErr = TMath::Sqrt(con*(1.-con)/kaonCon[i+1]);
 493     effKa->SetPoint(i,pt[i],(Float_t)kaonEff[i+1]/TMath::Max(kaonTot[i+1],1));
 494     effKa->SetPointError(i,0,effErr);
 495     conKa->SetPoint(i,pt[i],(Float_t)kaonNot[i+1]/TMath::Max(kaonCon[i+1],1));
 496     conKa->SetPointError(i,0,conErr);
 497   }
 498
 499   TCanvas *pCan=new TCanvas("Hmp","Efficiency and contamination",500,900);
 500   pCan->Divide(1,3);
 501
 502   pCan->cd(1);
 503   effPi->SetTitle("Pions");
 504   effPi->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 505   effPi->GetYaxis()->SetRangeUser(0.,1.);
 506   effPi->SetMarkerStyle(20);
 507   effPi->Draw("ALP");
 508   conPi->SetMarkerStyle(21);
 509   conPi->Draw("sameLP");
 510
 511   pCan->cd(2);
 512   effKa->SetTitle("Kaons");
 513   effKa->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 514   effKa->GetYaxis()->SetRangeUser(0.,1.);
 515   effKa->SetMarkerStyle(20);
 516   effKa->Draw("ALP");
 517   conKa->SetMarkerStyle(21);
 518   conKa->Draw("sameLP");
 519
 520   pCan->cd(3);
 521   effPr->SetTitle("Protons");
 522   effPr->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 523   effPr->GetYaxis()->SetRangeUser(0.,1.);
 524   effPr->SetMarkerStyle(20);
 525   effPr->Draw("ALP");
 526   conPr->SetMarkerStyle(21);
 527   conPr->Draw("sameLP");
 528
 529   TFile *outFile = new TFile("HmpidGraphs.root","recreate");
 530   pCan->Write();
 531   outFile->Close();
 532
 533   AliAnalysisTaskSE::Terminate();
 534
 535 }
 536
 537 //___________________________________________________________________________
 538 void AliHMPIDAnalysisTask::UserCreateOutputObjects() {
 539   //
 540   //HERE ONE CAN CREATE OUTPUT OBJECTS
 541   //TO BE SET BEFORE THE EXECUTION OF THE TASK
 542   //
 543
 544   //slot #1
 545    OpenFile(1);
 546    fHmpHistList = new TList();
 547
 548    fHmpPesdPhmp = new TH2F("fHmpPesdPhmp","HMPID: ESD p - running p;HMP p (GeV/c);ESD p (Gev/c)",100,0,10,100,0,10);
 549    fHmpHistList->Add(fHmpPesdPhmp);
 550
 551    fHmpCkovPesd = new TH2F("fHmpCkovPesd","HMPID: ThetaCherenkov vs P;p_esd (GeV/c);#Theta_C;Entries",100,0,10,110,0,1.1);
 552    fHmpHistList->Add(fHmpCkovPesd);
 553
 554    fHmpCkovPhmp = new TH2F("fHmpCkovPhmp","HMPID: ThetaCherenkov vs P;p_hmp (GeV/c);#Theta_C;Entries",100,0,10,110,0,1.1);
 555    fHmpHistList->Add(fHmpCkovPhmp);
 556
 557    fHmpMipTrkDist = new TH1F("fHmpMipTrkDist","HMPID MIP-Track distance;distance (cm);Entries",800,-20,20);
 558    fHmpHistList->Add(fHmpMipTrkDist);
 559
 560    fHmpMipTrkDistX = new TH1F("fHmpMipTrkDistX","HMPID MIP-Track distance in local X;distance (cm);Entries",800,-20,20);
 561    fHmpHistList->Add(fHmpMipTrkDistX);
 562
 563    fHmpMipTrkDistY = new TH1F("fHmpMipTrkDistY","HMPID MIP-Track distance in local Y;distance (cm);Entries",800,-20,20);
 564    fHmpHistList->Add(fHmpMipTrkDistY);
 565
 566    fHmpMipCharge3cm = new TH1F("fHmpMipCharge3cm","HMPID MIP Charge;MIP Charge (ADC);Entries",5001,-0.5,5000.5);
 567    fHmpHistList->Add(fHmpMipCharge3cm);
 568
 569    fHmpMipCharge1cm = new TH1F("fHmpMipCharge1cm","HMPID MIP Charge;MIP Charge (ADC);Entries",5001,-0.5,5000.5);
 570    fHmpHistList->Add(fHmpMipCharge1cm);
 571
 572    fHmpNumPhots = new TH1F("fHmpNumPhots","HMPID Number of photon clusters on ring;#photon clus.;Entries",51,-0.5,50.5);
 573    fHmpHistList->Add(fHmpNumPhots);
 574
 575    fHmpTrkFlags = new TH1F("fHmpTrkFlags","HMPID track flags",6,0,6);
 576    TString summary[6] =  {"NotPerformed","MipDistCut", "MipQdcCut", "NoPhotAccept", "kNoRad", "other"};
 577    for(Int_t ibin = 0; ibin < 6; ibin++) fHmpTrkFlags->GetXaxis()->SetBinLabel(ibin+1,Form("%i  %s",ibin+1,summary[ibin].Data()));
 578    fHmpHistList->Add(fHmpTrkFlags);
 579
 580    fPionEff = new TH1F("PionEff","Identified pions",fN1,0,fN1);
 581    fKaonEff = new TH1F("KaonEff","Identified kaons",fN1,0,fN1);
 582    fProtEff = new TH1F("ProtEff","Identified protons",fN2,0,fN2);
 583    fPionTot = new TH1I("PionTot","Total MC pions",fN1,0,fN1);
 584    fKaonTot = new TH1I("KaonTot","Total MC kaons",fN1,0,fN1);
 585    fProtTot = new TH1I("ProtTot","Total MC protons",fN2,0,fN2);
 586    fPionNot = new TH1F("PionNot","Misidentified pions",fN1,0,fN1);
 587    fKaonNot = new TH1F("KaonNot","Misidentified kaons",fN1,0,fN1);
 588    fProtNot = new TH1F("ProtNot","Misidentified protons",fN2,0,fN2);
 589    fPionCon = new TH1I("PionCon","Total not MC pions",fN1,0,fN1);
 590    fKaonCon = new TH1I("KaonCon","Total not MC kaons",fN1,0,fN1);
 591    fProtCon = new TH1I("ProtCon","Total not MC protons",fN2,0,fN2);
 592
 593    fHmpHistList->Add(fPionEff); fHmpHistList->Add(fKaonEff); fHmpHistList->Add(fProtEff);
 594    fHmpHistList->Add(fPionTot); fHmpHistList->Add(fKaonTot); fHmpHistList->Add(fProtTot);
 595    fHmpHistList->Add(fPionNot); fHmpHistList->Add(fKaonNot); fHmpHistList->Add(fProtNot);
 596    fHmpHistList->Add(fPionCon); fHmpHistList->Add(fKaonCon); fHmpHistList->Add(fProtCon);
 597
 598    fThetavsPiFromK= new TH2F("ThetavsPiFromK","Theta vs p of pions from K;p_esd (GeV/c);#Theta_C",100,0,10,110,0,1.1);
 599    fHmpHistList->Add(fThetavsPiFromK);
 600
 601    fThetapivsPesd = new TH2F("ThetapivsPesd","Theta of pions vs p of esd;p_esd (GeV/c);#Theta_C",100,0,10,110,0,1.1);
 602    fHmpHistList->Add(fThetapivsPesd);
 603
 604    fThetaKvsPesd  = new TH2F("ThetaKvsPesd","Theta of kaons vs p of esd;p_esd (GeV/c);#Theta_C",100,0,10,110,0,1.1);
 605    fHmpHistList->Add(fThetaKvsPesd);
 606
 607    fThetaPvsPesd  = new TH2F("ThetaPvsPesd","Theta of protons vs p of esd;p_esd (GeV/c);#Theta_C",100,0,10,110,0,1.1);
 608    fHmpHistList->Add(fThetaPvsPesd);
 609
 610    OpenFile(2);
 611    fTree = new TTree("Tree","Tree with data");
 612    fTree->Branch("Chamber",&fVar[0]);
 613    fTree->Branch("pHmp3",&fVar[1]);
 614    fTree->Branch("P",&fVar[2]);
 615    fTree->Branch("Xpc",&fVar[3]);
 616    fTree->Branch("Ypc",&fVar[4]);
 617    fTree->Branch("X",&fVar[5]);
 618    fTree->Branch("Y",&fVar[6]);
 619    fTree->Branch("ThetaTheorPion",&fVar[7]);
 620    fTree->Branch("ThetaTheorKaon",&fVar[8]);
 621    fTree->Branch("ThetaTheorProt",&fVar[9]);
 622    fTree->Branch("ThetaHmpTheorPion",&fVar[10]);
 623    fTree->Branch("ThetaHmpTheorKaon",&fVar[11]);
 624    fTree->Branch("ThetaHmpTheorProt",&fVar[12]);
 625    fTree->Branch("HMPIDsignal",&fVar[13]);
 626    fTree->Branch("Charge",&fVar[14]);
 627    fTree->Branch("Theta",&fVar[15]);
 628    fTree->Branch("Phi",&fVar[16]);
 629    fTree->Branch("Sign",&fVar[17]);
 630    fTree->Branch("NumPhotons",&fVar[18]);
 631    fTree->Branch("NumTPCclust",&fVar[19]);
 632    fTree->Branch("Prob0",&fVar[20]);
 633    fTree->Branch("Prob1",&fVar[21]);
 634    fTree->Branch("Prob2",&fVar[22]);
 635    fTree->Branch("Prob3",&fVar[23]);
 636    fTree->Branch("Prob4",&fVar[24]);
 637    fTree->Branch("TOFsignal",&fVar[25]);
 638    fTree->Branch("KinkIndex",&fVar[26]);
 639    fTree->Branch("Xv",&fVar[27]);
 640    fTree->Branch("Yv",&fVar[28]);
 641    fTree->Branch("Zv",&fVar[29]);
 642    fTree->Branch("TPCchi2",&fVar[30]);
 643    fTree->Branch("b0",&fVar[31]);
 644    fTree->Branch("b1",&fVar[32]);
 645    fTree->Branch("ClustSize",&fVar[33]);
 646 }
 647
 648 //____________________________________________________________________________________________________________________________________
 649 Bool_t AliHMPIDAnalysisTask::Equal(Double_t x, Double_t y, Double_t tolerance)
 650 {
 651  return abs(x - y) <= tolerance ;
 652 }
 653
 654 #endif