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