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       Double_t thetaCkov = track->GetHMPIDsignal() - (Int_t)track->GetHMPIDsignal();
 315       if (pHmp3) fHmpPesdPhmp->Fill(track->P(),pHmp3);
 316       if (dist<=1.0) fHmpMipCharge1cm->Fill(q);
 317       fHmpNumPhots->Fill(nph);
 318       fHmpCkovPesd->Fill(track->P(),thetaCkov);
 319       if (pHmp3) fHmpCkovPhmp->Fill(pHmp3,thetaCkov);
 320
 321       if (fUseMC && dist<distCut && TMath::Abs(th)<thetaCut){
 322         if (!pStack->IsPhysicalPrimary(label)) continue;
 323         Int_t pdgCode = TMath::Abs(pPart->GetPdgCode());
 324         if (pdgCode==211){
 325           fThetapivsPesd->Fill(track->P(),thetaCkov);
 326           Int_t mot=pPart->GetFirstMother();
 327           if (mot > -1){
 328             TParticle *pMot=pStack->Particle(mot);
 329             TString str=pMot->GetName();
 330             if (str.Contains("K")) fThetavsPiFromK->Fill(pHmp3,thetaCkov);
 331           }
 332         }
 333         if (pdgCode==321) fThetaKvsPesd->Fill(track->P(),thetaCkov);
 334         if (pdgCode==2212) fThetaPvsPesd->Fill(track->P(),thetaCkov);
 335
 336         if (track->P()<1. || track->P()>5.) continue;
 337         Int_t pBin=(Int_t) (2*(track->P()-1));
 338         if (pdgCode!=2212) fProtCon->Fill(pBin);
 339         if (pdgCode==2212){
 340           fProtTot->Fill(pBin);
 341           fProtEff->Fill(pBin,pPid->GetProbability(AliPID::kProton));
 342           fPionNot->Fill(pBin,pPid->GetProbability(AliPID::kPion));
 343           fKaonNot->Fill(pBin,pPid->GetProbability(AliPID::kKaon));
 344         }
 345         if (pdgCode!=211) fPionCon->Fill(pBin);
 346         if (pdgCode!=321) fKaonCon->Fill(pBin);
 347         if (pdgCode==211){
 348           if (pBin < 6){
 349             Float_t weight=pPid->GetProbability(AliPID::kElectron)+
 350                            pPid->GetProbability(AliPID::kMuon)+
 351                            pPid->GetProbability(AliPID::kPion);
 352             fPionTot->Fill(pBin);
 353             fPionEff->Fill(pBin,weight);
 354             fKaonNot->Fill(pBin,pPid->GetProbability(AliPID::kKaon));
 355           }
 356           fProtNot->Fill(pBin,pPid->GetProbability(AliPID::kProton));
 357         }
 358         if (pdgCode==321){
 359           if (pBin < 6){
 360             fKaonTot->Fill(pBin);
 361             fKaonEff->Fill(pBin,pPid->GetProbability(AliPID::kKaon));
 362             fPionNot->Fill(pBin,(pPid->GetProbability(AliPID::kPion)));
 363           }
 364           fProtNot->Fill(pBin,(pPid->GetProbability(AliPID::kProton)));
 365         }
 366       }
 367     }//there is signal
 368     fVar[0] = track->GetHMPIDcluIdx()/1000000;
 369     fVar[1] = pHmp3;
 370     fVar[2] = (Float_t)track->P();
 371     fVar[3] = xpc;
 372     fVar[4] = ypc;
 373     fVar[5] = x;
 374     fVar[6] = y;
 375     fVar[7] = (Float_t)thetaCkov;
 376     fVar[8] = q;
 377     fVar[9] = th;
 378     fVar[10] = ph;
 379     fVar[11] = (Float_t)track->GetSign();
 380     fVar[12] = (Float_t)nph;
 381     fVar[13] = (Float_t)track->GetNcls(1);
 382     fVar[14] = (Float_t)probs[0];
 383     fVar[15] = (Float_t)probs[1];
 384     fVar[16] = (Float_t)probs[2];
 385     fVar[17] = (Float_t)probs[3];
 386     fVar[18] = (Float_t)probs[4];
 387     fVar[19] = (Float_t)track->GetTOFsignal();
 388     fVar[20] = (Float_t)track->GetKinkIndex(0);
 389     fVar[21] = (Float_t)track->Xv();
 390     fVar[22] = (Float_t)track->Yv();
 391     fVar[23] = (Float_t)track->Zv();
 392     fVar[24] = (Float_t)track->GetTPCchi2();
 393     fVar[25] = b[0];
 394     fVar[26] = b[1];
 395     fVar[27] = track->GetHMPIDcluIdx()%1000000/1000;
 396     fTree->Fill();
 397   }//track loop
 398   delete pPid;
 399
 400   if (fUseMC){
 401     Float_t phiMin=(5./180.)*TMath::Pi() , phiMax=(55./180.)*TMath::Pi();
 402     for (Int_t iPart=0; iPart<pStack->GetNprimary(); iPart++){
 403       TString namepart=pStack->Particle(iPart)->GetName();
 404       if (!namepart.Contains("proton")) continue;
 405       pPart=pStack->Particle(iPart);
 406       if (!pStack->IsPhysicalPrimary(iPart) || pPart->P()<1. || pPart->P()>5.) continue;
 407       if (pPart->Phi()< phiMin || pPart->Phi()>phiMax) continue;
 408       if (TMath::Abs(pPart->Eta()) > 0.5) continue;
 409       Int_t pBin=(Int_t) (2*(pPart->P()-1));
 410       if (namepart == "proton") fProtGen->Fill(pBin);
 411       if (namepart == "antiproton") fPbarGen->Fill(pBin);
 412       for (Int_t iTrack=0; iTrack<fESD->GetNumberOfTracks(); iTrack++){
 413         if (fESD->GetTrack(iTrack)->GetLabel() != iPart) continue;
 414         if (fESD->GetTrack(iTrack)->GetHMPIDsignal() < 0 ) continue;
 415         if (namepart == "proton") fProtHmp->Fill(pBin);
 416         if (namepart == "antiproton") fPbarHmp->Fill(pBin);
 417       }
 418     }
 419   }
 420
 421   /* PostData(0) is taken care of by AliAnalysisTaskSE */
 422   PostData(1,fHmpHistList);
 423   PostData(2,fTree);
 424 }
 425
 426 //___________________________________________________________________________
 427 void AliHMPIDAnalysisTask::Terminate(Option_t*)
 428 {
 429   // The Terminate() function is the last function to be called during
 430   // a query. It always runs on the client, it can be used to present
 431   // the results graphically or save the results to file.
 432
 433   Info("Terminate"," ");
 434
 435   if (!fUseMC) return;
 436
 437   fHmpHistList = dynamic_cast<TList*> (GetOutputData(1));
 438
 439   if (!fHmpHistList) {
 440     AliError("Histogram List is not available");
 441     return;
 442   }
 443
 444   fPionEff = dynamic_cast<TH1F*> (fHmpHistList->FindObject("PionEff"));
 445   fPionTot = dynamic_cast<TH1I*> (fHmpHistList->FindObject("PionTot"));
 446   fPionNot = dynamic_cast<TH1F*> (fHmpHistList->FindObject("PionNot"));
 447   fPionCon = dynamic_cast<TH1I*> (fHmpHistList->FindObject("PionCon"));
 448   fKaonEff = dynamic_cast<TH1F*> (fHmpHistList->FindObject("KaonEff"));
 449   fKaonTot = dynamic_cast<TH1I*> (fHmpHistList->FindObject("KaonTot"));
 450   fKaonNot = dynamic_cast<TH1F*> (fHmpHistList->FindObject("KaonNot"));
 451   fKaonCon = dynamic_cast<TH1I*> (fHmpHistList->FindObject("KaonCon"));
 452   fProtEff = dynamic_cast<TH1F*> (fHmpHistList->FindObject("ProtEff"));
 453   fProtTot = dynamic_cast<TH1I*> (fHmpHistList->FindObject("ProtTot"));
 454   fProtNot = dynamic_cast<TH1F*> (fHmpHistList->FindObject("ProtNot"));
 455   fProtCon = dynamic_cast<TH1I*> (fHmpHistList->FindObject("ProtCon"));
 456
 457   Float_t *pionEff = fPionEff->GetArray();
 458   Int_t   *pionTot = fPionTot->GetArray();
 459   Float_t *pionNot = fPionNot->GetArray();
 460   Int_t   *pionCon = fPionCon->GetArray();
 461   Float_t *kaonEff = fKaonEff->GetArray();
 462   Int_t   *kaonTot = fKaonTot->GetArray();
 463   Float_t *kaonNot = fKaonNot->GetArray();
 464   Int_t   *kaonCon = fKaonCon->GetArray();
 465   Float_t *protEff = fProtEff->GetArray();
 466   Int_t   *protTot = fProtTot->GetArray();
 467   Float_t *protNot = fProtNot->GetArray();
 468   Int_t   *protCon = fProtCon->GetArray();
 469
 470   TGraphErrors *effPi = new TGraphErrors(fN1);
 471   TGraphErrors *effKa = new TGraphErrors(fN1);
 472   TGraphErrors *effPr = new TGraphErrors(fN2);
 473   TGraphErrors *conPi = new TGraphErrors(fN1);
 474   TGraphErrors *conKa = new TGraphErrors(fN1);
 475   TGraphErrors *conPr = new TGraphErrors(fN2);
 476
 477   Float_t pt[8]={1.25,1.75,2.25,2.75,3.25,3.75,4.25,4.75};
 478   Float_t eff=0, effErr=0, con=0, conErr=0;
 479   for (Int_t i=0; i< fN2; i++){
 480     eff = protEff[i+1]/TMath::Max(protTot[i+1],1);
 481     effErr = TMath::Sqrt(eff*(1.-eff)/TMath::Max(protTot[i+1],1));
 482     con = protNot[i+1]/TMath::Max(protCon[i+1],1);
 483     conErr = TMath::Sqrt(con*(1.-con)/protCon[i+1]);
 484     effPr->SetPoint(i,pt[i],eff);
 485     effPr->SetPointError(i,0,effErr);
 486     conPr->SetPoint(i,pt[i],con);
 487     conPr->SetPointError(i,0,conErr);
 488
 489     if (i>=fN1) continue;
 490     eff = pionEff[i+1]/pionTot[i+1];
 491     effErr = TMath::Sqrt(eff*(1.-eff)/pionTot[i+1]);
 492     con=pionNot[i+1]/pionCon[i+1];
 493     conErr = TMath::Sqrt(con*(1.-con)/pionCon[i+1]);
 494     effPi->SetPoint(i,pt[i],(Float_t)pionEff[i+1]/(Float_t)pionTot[i+1]);
 495     effPi->SetPointError(i,0,effErr);
 496     conPi->SetPoint(i,pt[i],(Float_t)pionNot[i+1]/(Float_t)pionCon[i+1]);
 497     conPi->SetPointError(i,0,conErr);
 498
 499     eff = kaonEff[i+1]/TMath::Max(kaonTot[i+1],1);
 500     effErr = TMath::Sqrt(eff*(1.-eff)/kaonTot[i+1]);
 501     con = kaonNot[i+1]/TMath::Max(kaonCon[i+1],1);
 502     conErr = TMath::Sqrt(con*(1.-con)/kaonCon[i+1]);
 503     effKa->SetPoint(i,pt[i],(Float_t)kaonEff[i+1]/TMath::Max(kaonTot[i+1],1));
 504     effKa->SetPointError(i,0,effErr);
 505     conKa->SetPoint(i,pt[i],(Float_t)kaonNot[i+1]/TMath::Max(kaonCon[i+1],1));
 506     conKa->SetPointError(i,0,conErr);
 507   }
 508
 509   fProtGen = dynamic_cast<TH1I*> (fHmpHistList->FindObject("ProtGen"));
 510   fPbarGen = dynamic_cast<TH1I*> (fHmpHistList->FindObject("PbarGen"));
 511   fProtHmp = dynamic_cast<TH1I*> (fHmpHistList->FindObject("ProtHmp"));
 512   fPbarHmp = dynamic_cast<TH1I*> (fHmpHistList->FindObject("PbarHmp"));
 513
 514   Int_t   *protGen = fProtGen->GetArray();
 515   Int_t   *pbarGen = fPbarGen->GetArray();
 516   Int_t   *protHmp = fProtHmp->GetArray();
 517   Int_t   *pbarHmp = fPbarHmp->GetArray();
 518
 519   TGraphErrors *ratPr = new TGraphErrors(fN2);
 520   TGraphErrors *ratPb = new TGraphErrors(fN2);
 521
 522   Float_t rat1=0, rat1Err=0, rat2=0, rat2Err=0;
 523   for (Int_t i=0; i< fN2; i++){
 524     rat1 = (Float_t)protHmp[i+1]/TMath::Max(protGen[i+1],1);
 525     rat1Err = TMath::Sqrt(rat1*(1.-rat1)/TMath::Max(protGen[i+1],1));
 526     rat2 = (Float_t)pbarHmp[i+1]/TMath::Max(pbarGen[i+1],1);
 527     rat2Err = TMath::Sqrt(rat2*(1.-rat2)/TMath::Max(pbarGen[i+1],1));
 528     ratPr->SetPoint(i,pt[i],rat1);
 529     ratPr->SetPointError(i,0,rat1Err);
 530     ratPb->SetPoint(i,pt[i],rat2);
 531     ratPb->SetPointError(i,0,rat2Err);
 532   }
 533
 534   TCanvas *pCan=new TCanvas("Hmp","Efficiency and contamination",500,900);
 535   pCan->Divide(2,3);
 536
 537   pCan->cd(1);
 538   effPi->SetTitle("Pions");
 539   effPi->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 540   effPi->GetYaxis()->SetRangeUser(0.,1.);
 541   effPi->SetMarkerStyle(20);
 542   effPi->Draw("ALP");
 543   conPi->SetMarkerStyle(21);
 544   conPi->Draw("sameLP");
 545
 546   pCan->cd(3);
 547   effKa->SetTitle("Kaons");
 548   effKa->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 549   effKa->GetYaxis()->SetRangeUser(0.,1.);
 550   effKa->SetMarkerStyle(20);
 551   effKa->Draw("ALP");
 552   conKa->SetMarkerStyle(21);
 553   conKa->Draw("sameLP");
 554
 555   pCan->cd(5);
 556   effPr->SetTitle("Protons");
 557   effPr->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 558   effPr->GetYaxis()->SetRangeUser(0.,1.);
 559   effPr->SetMarkerStyle(20);
 560   effPr->Draw("ALP");
 561   conPr->SetMarkerStyle(21);
 562   conPr->Draw("sameLP");
 563
 564   pCan->cd(2);
 565   ratPr->SetTitle("Ratios");
 566   ratPr->GetXaxis()->SetTitle("p_{T} (GeV/c)");
 567   ratPr->GetYaxis()->SetRangeUser(0.,1.);
 568   ratPr->SetMarkerStyle(20);
 569   ratPr->Draw("ALP");
 570   ratPb->SetMarkerStyle(21);
 571   ratPb->Draw("sameLP");
 572
 573   TFile *outFile = new TFile("HmpidGraphs.root","recreate");
 574   pCan->Write();
 575   outFile->Close();
 576
 577   AliAnalysisTaskSE::Terminate();
 578
 579 }
 580
 581 //___________________________________________________________________________
 582 void AliHMPIDAnalysisTask::UserCreateOutputObjects() {
 583   //
 584   //HERE ONE CAN CREATE OUTPUT OBJECTS
 585   //TO BE SET BEFORE THE EXECUTION OF THE TASK
 586   //
 587
 588   //slot #1
 589 //   OpenFile(1);
 590    fHmpHistList = new TList();
 591    fHmpHistList->SetOwner();
 592
 593    fHmpPesdPhmp = new TH2F("fHmpPesdPhmp","HMPID: ESD p - running p;HMP p (GeV/c);ESD p (Gev/c)",100,0,10,100,0,10);
 594    fHmpHistList->Add(fHmpPesdPhmp);
 595
 596    fHmpCkovPesd = new TH2F("fHmpCkovPesd","HMPID: ThetaCherenkov vs P;p_esd (GeV/c);#Theta_C;Entries",100,0,10,110,0,1.1);
 597    fHmpHistList->Add(fHmpCkovPesd);
 598
 599    fHmpCkovPhmp = new TH2F("fHmpCkovPhmp","HMPID: ThetaCherenkov vs P;p_hmp (GeV/c);#Theta_C;Entries",100,0,10,110,0,1.1);
 600    fHmpHistList->Add(fHmpCkovPhmp);
 601
 602    fHmpMipTrkDist = new TH1F("fHmpMipTrkDist","HMPID MIP-Track distance;distance (cm);Entries",400,0,20);
 603    fHmpHistList->Add(fHmpMipTrkDist);
 604
 605    fHmpMipTrkDistX = new TH1F("fHmpMipTrkDistX","HMPID MIP-Track distance in local X;distance (cm);Entries",800,-20,20);
 606    fHmpHistList->Add(fHmpMipTrkDistX);
 607
 608    fHmpMipTrkDistY = new TH1F("fHmpMipTrkDistY","HMPID MIP-Track distance in local Y;distance (cm);Entries",800,-20,20);
 609    fHmpHistList->Add(fHmpMipTrkDistY);
 610
 611    fHmpMipCharge3cm = new TH1F("fHmpMipCharge3cm","HMPID MIP Charge;MIP Charge (ADC);Entries",5001,-0.5,5000.5);
 612    fHmpHistList->Add(fHmpMipCharge3cm);
 613
 614    fHmpMipCharge1cm = new TH1F("fHmpMipCharge1cm","HMPID MIP Charge;MIP Charge (ADC);Entries",5001,-0.5,5000.5);
 615    fHmpHistList->Add(fHmpMipCharge1cm);
 616
 617    fHmpNumPhots = new TH1F("fHmpNumPhots","HMPID Number of photon clusters on ring;#photon clus.;Entries",51,-0.5,50.5);
 618    fHmpHistList->Add(fHmpNumPhots);
 619
 620    fHmpTrkFlags = new TH1F("fHmpTrkFlags","HMPID track flags",6,0,6);
 621    TString summary[6] =  {"NotPerformed","MipDistCut", "MipQdcCut", "NoPhotAccept", "kNoRad", "other"};
 622    for(Int_t ibin = 0; ibin < 6; ibin++) fHmpTrkFlags->GetXaxis()->SetBinLabel(ibin+1,Form("%i  %s",ibin+1,summary[ibin].Data()));
 623    fHmpHistList->Add(fHmpTrkFlags);
 624
 625    fPionEff = new TH1F("PionEff","Identified pions",fN1,0,fN1);
 626    fKaonEff = new TH1F("KaonEff","Identified kaons",fN1,0,fN1);
 627    fProtEff = new TH1F("ProtEff","Identified protons",fN2,0,fN2);
 628    fPionTot = new TH1I("PionTot","Total MC pions",fN1,0,fN1);
 629    fKaonTot = new TH1I("KaonTot","Total MC kaons",fN1,0,fN1);
 630    fProtTot = new TH1I("ProtTot","Total MC protons",fN2,0,fN2);
 631    fPionNot = new TH1F("PionNot","Misidentified pions",fN1,0,fN1);
 632    fKaonNot = new TH1F("KaonNot","Misidentified kaons",fN1,0,fN1);
 633    fProtNot = new TH1F("ProtNot","Misidentified protons",fN2,0,fN2);
 634    fPionCon = new TH1I("PionCon","Total not MC pions",fN1,0,fN1);
 635    fKaonCon = new TH1I("KaonCon","Total not MC kaons",fN1,0,fN1);
 636    fProtCon = new TH1I("ProtCon","Total not MC protons",fN2,0,fN2);
 637
 638    fHmpHistList->Add(fPionEff); fHmpHistList->Add(fKaonEff); fHmpHistList->Add(fProtEff);
 639    fHmpHistList->Add(fPionTot); fHmpHistList->Add(fKaonTot); fHmpHistList->Add(fProtTot);
 640    fHmpHistList->Add(fPionNot); fHmpHistList->Add(fKaonNot); fHmpHistList->Add(fProtNot);
 641    fHmpHistList->Add(fPionCon); fHmpHistList->Add(fKaonCon); fHmpHistList->Add(fProtCon);
 642
 643    fProtGen = new TH1I("ProtGen","Generated protons",fN2,0,fN2);
 644    fPbarGen = new TH1I("PbarGen","Generated antiprotons",fN2,0,fN2);
 645    fProtHmp = new TH1I("ProtHmp","HMPID protons",fN2,0,fN2);
 646    fPbarHmp = new TH1I("PbarHmp","HMPID antiprotons",fN2,0,fN2);
 647
 648    fHmpHistList->Add(fProtGen); fHmpHistList->Add(fPbarGen);
 649    fHmpHistList->Add(fProtHmp); fHmpHistList->Add(fPbarHmp);
 650
 651    fThetavsPiFromK= new TH2F("ThetavsPiFromK","Theta vs p of pions from K;p_esd (GeV/c);#Theta_C",100,0,10,110,0,1.1);
 652    fHmpHistList->Add(fThetavsPiFromK);
 653
 654    fThetapivsPesd = new TH2F("ThetapivsPesd","Theta of pions vs p of esd;p_esd (GeV/c);#Theta_C",100,0,10,110,0,1.1);
 655    fHmpHistList->Add(fThetapivsPesd);
 656
 657    fThetaKvsPesd  = new TH2F("ThetaKvsPesd","Theta of kaons vs p of esd;p_esd (GeV/c);#Theta_C",100,0,10,110,0,1.1);
 658    fHmpHistList->Add(fThetaKvsPesd);
 659
 660    fThetaPvsPesd  = new TH2F("ThetaPvsPesd","Theta of protons vs p of esd;p_esd (GeV/c);#Theta_C",100,0,10,110,0,1.1);
 661    fHmpHistList->Add(fThetaPvsPesd);
 662
 663 //   OpenFile(2);
 664    fTree = new TTree("Tree","Tree with data");
 665    fTree->Branch("Chamber",&fVar[0]);
 666    fTree->Branch("pHmp3",&fVar[1]);
 667    fTree->Branch("P",&fVar[2]);
 668    fTree->Branch("Xpc",&fVar[3]);
 669    fTree->Branch("Ypc",&fVar[4]);
 670    fTree->Branch("X",&fVar[5]);
 671    fTree->Branch("Y",&fVar[6]);
 672    fTree->Branch("HMPIDsignal",&fVar[7]);
 673    fTree->Branch("Charge",&fVar[8]);
 674    fTree->Branch("Theta",&fVar[9]);
 675    fTree->Branch("Phi",&fVar[10]);
 676    fTree->Branch("Sign",&fVar[11]);
 677    fTree->Branch("NumPhotons",&fVar[12]);
 678    fTree->Branch("NumTPCclust",&fVar[13]);
 679    fTree->Branch("Prob0",&fVar[14]);
 680    fTree->Branch("Prob1",&fVar[15]);
 681    fTree->Branch("Prob2",&fVar[16]);
 682    fTree->Branch("Prob3",&fVar[17]);
 683    fTree->Branch("Prob4",&fVar[18]);
 684    fTree->Branch("TOFsignal",&fVar[19]);
 685    fTree->Branch("KinkIndex",&fVar[20]);
 686    fTree->Branch("Xv",&fVar[21]);
 687    fTree->Branch("Yv",&fVar[22]);
 688    fTree->Branch("Zv",&fVar[23]);
 689    fTree->Branch("TPCchi2",&fVar[24]);
 690    fTree->Branch("b0",&fVar[25]);
 691    fTree->Branch("b1",&fVar[26]);
 692    fTree->Branch("ClustSize",&fVar[27]);
 693
 694    PostData(1,fHmpHistList);
 695    PostData(2,fTree);
 696 }
 697
 698 //____________________________________________________________________________________________________________________________________
 699 Bool_t AliHMPIDAnalysisTask::Equal(Double_t x, Double_t y, Double_t tolerance)
 700 {
 701  return abs(x - y) <= tolerance ;
 702 }
 703
 704 #endif