PWG1/HMPID/AliHMPIDTaskQA.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 // AliHMPIDTaskQA - Class representing a quality check tool of HMPID
  18 // A set of histograms is created.
  19 //==============================================================================
  20
  21 #ifndef AliHMPIDTaskQA_CXX
  22 #define AliHMPIDTaskQA_CXX
  23
  24
  25 #include "TH1.h"
  26 #include "TH2.h"
  27 #include "TFile.h"
  28 #include "TCanvas.h"
  29 #include "TGraphErrors.h"
  30 #include "AliAnalysisManager.h"
  31 #include "AliESDInputHandler.h"
  32 #include "AliMCEventHandler.h"
  33 #include "AliMCEvent.h"
  34 #include "AliESDtrack.h"
  35 #include "AliPID.h"
  36 #include "AliLog.h"
  37 #include "AliHMPIDTaskQA.h"
  38
  39 ClassImp(AliHMPIDTaskQA)
  40
  41 //__________________________________________________________________________
  42 AliHMPIDTaskQA::AliHMPIDTaskQA() :
  43   fESD(0x0),fMC(0x0),fUseMC(kTRUE),
  44   fHmpHistList(0x0),
  45   fHmpPesdPhmp(0x0),fHmpCkovPesd(0x0),fHmpCkovPhmp(0x0),
  46   fHmpMipTrkDistX(0x0),fHmpMipCharge3cm(0x0),
  47   fHmpTrkFlags(0x0),
  48   fN1(6),
  49   fN2(8),
  50   fPionEff(0x0),
  51   fKaonEff(0x0),
  52   fProtEff(0x0),
  53   fPionTot(0x0),
  54   fKaonTot(0x0),
  55   fProtTot(0x0),
  56   fPionNot(0x0),
  57   fKaonNot(0x0),
  58   fProtNot(0x0),
  59   fPionCon(0x0),
  60   fKaonCon(0x0),
  61   fProtCon(0x0)
  62 {
  63   //
  64   //Default ctor
  65   //
  66   for (Int_t i=0; i<7; i++)
  67     fHmpMipTrkDistPosY[i] = fHmpMipTrkDistNegY[i] = 0x0;
  68 }
  69
  70 //___________________________________________________________________________
  71 AliHMPIDTaskQA::AliHMPIDTaskQA(const Char_t* name) :
  72   AliAnalysisTaskSE(name),
  73   fESD(0x0), fMC(0x0), fUseMC(kTRUE),
  74   fHmpHistList(0x0),
  75   fHmpPesdPhmp(0x0),fHmpCkovPesd(0x0),fHmpCkovPhmp(0x0),
  76   fHmpMipTrkDistX(0x0),fHmpMipCharge3cm(0x0),
  77   fHmpTrkFlags(0x0),
  78   fN1(6),
  79   fN2(8),
  80   fPionEff(0x0),
  81   fKaonEff(0x0),
  82   fProtEff(0x0),
  83   fPionTot(0x0),
  84   fKaonTot(0x0),
  85   fProtTot(0x0),
  86   fPionNot(0x0),
  87   fKaonNot(0x0),
  88   fProtNot(0x0),
  89   fPionCon(0x0),
  90   fKaonCon(0x0),
  91   fProtCon(0x0)
  92 {
  93   //
  94   // Constructor. Initialization of Inputs and Outputs
  95   //
  96   for (Int_t i=0; i<7; i++)
  97     fHmpMipTrkDistPosY[i] = fHmpMipTrkDistNegY[i] = 0x0;
  98
  99   DefineOutput(1,TList::Class());
 100 }
 101
 102 //___________________________________________________________________________
 103 AliHMPIDTaskQA& AliHMPIDTaskQA::operator=(const AliHMPIDTaskQA& c)
 104 {
 105   //
 106   // Assignment operator
 107   //
 108   if (this!=&c) {
 109     AliAnalysisTaskSE::operator=(c);
 110     fESD             = c.fESD;
 111     fMC              = c.fMC;
 112     fUseMC           = c.fUseMC;
 113     fHmpHistList     = c.fHmpHistList;
 114     fHmpPesdPhmp     = c.fHmpPesdPhmp;
 115     fHmpCkovPesd     = c.fHmpCkovPesd;
 116     fHmpCkovPhmp     = c.fHmpCkovPhmp;
 117     fHmpMipTrkDistX  = c.fHmpMipTrkDistX;
 118     fHmpMipCharge3cm = c.fHmpMipCharge3cm;
 119     fHmpTrkFlags     = c.fHmpTrkFlags;
 120     fN1              = c.fN1;
 121     fN2              = c.fN2;
 122     fPionEff         = c.fPionEff;
 123     fKaonEff         = c.fKaonEff;
 124     fProtEff         = c.fProtEff;
 125     fPionTot         = c.fPionTot;
 126     fKaonTot         = c.fKaonTot;
 127     fProtTot         = c.fProtTot;
 128     fPionNot         = c.fPionNot;
 129     fKaonNot         = c.fKaonNot;
 130     fProtNot         = c.fProtNot;
 131     fPionCon         = c.fPionCon;
 132     fKaonCon         = c.fKaonCon;
 133     fProtCon         = c.fProtCon;
 134     for (Int_t i=0; i<7; i++){
 135       fHmpMipTrkDistPosY[i] = c.fHmpMipTrkDistPosY[i];
 136       fHmpMipTrkDistNegY[i] = c.fHmpMipTrkDistNegY[i];
 137     }
 138   }
 139   return *this;
 140 }
 141
 142 //___________________________________________________________________________
 143 AliHMPIDTaskQA::AliHMPIDTaskQA(const AliHMPIDTaskQA& c) :
 144   AliAnalysisTaskSE(c),
 145   fESD(c.fESD),fMC(c.fMC),fUseMC(c.fUseMC),
 146   fHmpHistList(c.fHmpHistList),
 147   fHmpPesdPhmp(c.fHmpPesdPhmp),fHmpCkovPesd(c.fHmpCkovPesd),fHmpCkovPhmp(c.fHmpCkovPhmp),
 148   fHmpMipTrkDistX(c.fHmpMipTrkDistX),
 149   fHmpMipCharge3cm(c.fHmpMipCharge3cm),
 150   fHmpTrkFlags(c.fHmpTrkFlags),
 151   fN1(c.fN1),
 152   fN2(c.fN2),
 153   fPionEff(c.fPionEff),
 154   fKaonEff(c.fKaonEff),
 155   fProtEff(c.fProtEff),
 156   fPionTot(c.fPionTot),
 157   fKaonTot(c.fKaonTot),
 158   fProtTot(c.fProtTot),
 159   fPionNot(c.fPionNot),
 160   fKaonNot(c.fKaonNot),
 161   fProtNot(c.fProtNot),
 162   fPionCon(c.fPionCon),
 163   fKaonCon(c.fKaonCon),
 164   fProtCon(c.fProtCon)
 165 {
 166   //
 167   // Copy Constructor
 168   //
 169   for (Int_t i=0; i<7; i++){
 170     fHmpMipTrkDistPosY[i] = c.fHmpMipTrkDistPosY[i];
 171     fHmpMipTrkDistNegY[i] = c.fHmpMipTrkDistNegY[i];
 172   }
 173 }
 174
 175 //___________________________________________________________________________
 176 AliHMPIDTaskQA::~AliHMPIDTaskQA() {
 177   //
 178   //destructor
 179   //
 180   Info("~AliHMPIDTaskQA","Calling Destructor");
 181   if (fHmpHistList /*&& !AliAnalysisManager::GetAnalysisManager()->IsProofMode()*/) delete fHmpHistList;
 182 }
 183
 184 //___________________________________________________________________________
 185 void AliHMPIDTaskQA::ConnectInputData(Option_t *option)
 186 {
 187   AliAnalysisTaskSE::ConnectInputData(option);
 188
 189   // Connect ESD here
 190   AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 191   if (!esdH) {
 192     AliDebug(2,Form("ERROR: Could not get ESDInputHandler"));
 193   } else
 194     fESD = esdH->GetEvent();
 195
 196   if (fUseMC){
 197     // Connect MC
 198     AliMCEventHandler *mcH = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 199     if (!mcH) {
 200       AliDebug(2,Form("ERROR: Could not get MCEventHandler"));
 201       fUseMC = kFALSE;
 202     } else
 203       fMC = mcH->MCEvent();
 204       if (!fMC) AliDebug(2,Form("ERROR: Could not get MCEvent"));
 205   }
 206 }
 207
 208 //___________________________________________________________________________
 209 void AliHMPIDTaskQA::UserExec(Option_t *)
 210 {
 211   Double_t priors[5]={1.,1.,1.,1.,1.}; //{0.01,0.01,0.83,0.10,0.5};
 212   Double_t probs[5];
 213   AliPID *pPid = new AliPID();
 214   pPid->SetPriors(priors);
 215   AliESDtrack *track=0;
 216   TParticle *pPart=0;
 217   AliStack* pStack = 0;
 218   Int_t label;
 219   if (fUseMC){
 220     pStack = fMC->Stack();
 221   }
 222
 223   //
 224   // Main loop function, executed on Event basis
 225   //
 226   for (Int_t iTrack = 0; iTrack<fESD->GetNumberOfTracks(); iTrack++) {
 227
 228     track = fESD->GetTrack(iTrack);
 229     if(!track) continue;
 230     Double_t rin[3], rout[3];
 231     track->GetInnerXYZ(rin);
 232     track->GetOuterXYZ(rout);
 233     Double_t ktol = 0.001;
 234
 235     if(Equal(track->GetHMPIDsignal(),-20.,ktol))      fHmpTrkFlags->Fill(0);
 236     else if(Equal(track->GetHMPIDsignal(),-9.,ktol))  fHmpTrkFlags->Fill(1);
 237     else if(Equal(track->GetHMPIDsignal(),-5.,ktol))  fHmpTrkFlags->Fill(2);
 238     else if(Equal(track->GetHMPIDsignal(),-11.,ktol)) fHmpTrkFlags->Fill(3);
 239     else if(Equal(track->GetHMPIDsignal(),-22.,ktol)) fHmpTrkFlags->Fill(4);
 240     else fHmpTrkFlags->Fill(4);
 241
 242     if(Equal(track->GetHMPIDsignal(),-20.,ktol)) continue;
 243     if(track->GetHMPIDcluIdx() < 0) continue;
 244
 245     Int_t q, nph;
 246     Float_t x, y;
 247     Float_t xpc, ypc, th, ph;
 248     track->GetHMPIDmip(x,y,q,nph);
 249     track->GetHMPIDtrk(xpc,ypc,th,ph);
 250
 251     if(Equal(x,0.,ktol) && Equal(y,0.,ktol) && Equal(xpc,0.,ktol) && Equal(ypc,0.,ktol)) continue;
 252
 253     Float_t sign = track->GetSign();
 254     Int_t ch = track->GetHMPIDcluIdx()/1000000;
 255     Double_t dist = TMath::Sqrt( (xpc-x)*(xpc-x) + (ypc - y)*(ypc - y));
 256     fHmpMipTrkDistX->Fill(xpc - x);
 257     if (sign > 0.) fHmpMipTrkDistPosY[ch]->Fill(ypc - y);
 258     if (sign < 0.) fHmpMipTrkDistNegY[ch]->Fill(ypc - y);
 259     Double_t pHmp[3] = {0}, pHmp3 = 0;
 260     if (track->GetOuterHmpPxPyPz(pHmp)) pHmp3 = TMath::Sqrt(pHmp[0]*pHmp[0]+pHmp[1]*pHmp[1]+pHmp[2]*pHmp[2]);
 261     if (dist <= 3.0) fHmpMipCharge3cm->Fill(q);
 262
 263     track->GetHMPIDpid(probs);
 264     pPid->SetProbabilities(probs);
 265     if (fUseMC){
 266       if ((label = track->GetLabel()) < 0) continue;
 267       pPart = pStack->Particle(label);
 268     }
 269
 270     if(track->GetHMPIDsignal() > 0 ){
 271       if (pHmp3) fHmpPesdPhmp->Fill(track->P(),pHmp3);
 272       fHmpCkovPesd->Fill(track->P(),track->GetHMPIDsignal());
 273       if (pHmp3) fHmpCkovPhmp->Fill(pHmp3,track->GetHMPIDsignal());
 274
 275       if (fUseMC && dist<0.5 && TMath::Abs(th)<0.13){
 276         if (!pStack->IsPhysicalPrimary(label)) continue;
 277         if (track->Pt()<1. || track->Pt()>5.) continue;
 278         Int_t pdgCode = TMath::Abs(pPart->GetPdgCode());
 279         Int_t ptBin=(Int_t) (2*(track->Pt()-1));
 280         if (pdgCode!=2212) fProtCon->Fill(ptBin);
 281         if (pdgCode==2212){
 282           fProtTot->Fill(ptBin);
 283           fProtEff->Fill(ptBin,pPid->GetProbability(AliPID::kProton));
 284           fPionNot->Fill(ptBin,pPid->GetProbability(AliPID::kPion));
 285           fKaonNot->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
 286         }
 287         if (pdgCode!=211) fPionCon->Fill(ptBin);
 288         if (pdgCode!=321) fKaonCon->Fill(ptBin);
 289         if (pdgCode==211){
 290           if (ptBin < 6){
 291             Float_t weight=pPid->GetProbability(AliPID::kElectron)+
 292                            pPid->GetProbability(AliPID::kMuon)+
 293                            pPid->GetProbability(AliPID::kPion);
 294             fPionTot->Fill(ptBin);
 295             fPionEff->Fill(ptBin,weight);
 296             fKaonNot->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
 297           }
 298           fProtNot->Fill(ptBin,pPid->GetProbability(AliPID::kProton));
 299         }
 300         if (pdgCode==321){
 301           if (ptBin < 6){
 302             fKaonTot->Fill(ptBin);
 303             fKaonEff->Fill(ptBin,pPid->GetProbability(AliPID::kKaon));
 304             fPionNot->Fill(ptBin,(pPid->GetProbability(AliPID::kPion)));
 305           }
 306           fProtNot->Fill(ptBin,(pPid->GetProbability(AliPID::kProton)));
 307         }
 308       }
 309     }//there is signal
 310   }//track loop
 311   delete pPid;
 312
 313   /* PostData(0) is taken care of by AliAnalysisTaskSE */
 314   PostData(1,fHmpHistList);
 315 }
 316
 317 //___________________________________________________________________________
 318 void AliHMPIDTaskQA::Terminate(Option_t*)
 319 {
 320   // The Terminate() function is the last function to be called during
 321   // a query. It always runs on the client, it can be used to present
 322   // the results graphically or save the results to file.
 323
 324   Info("Terminate"," ");
 325   AliAnalysisTaskSE::Terminate();
 326
 327 }
 328
 329 //___________________________________________________________________________
 330 void AliHMPIDTaskQA::UserCreateOutputObjects() {
 331   //
 332   //HERE ONE CAN CREATE OUTPUT OBJECTS
 333   //TO BE SET BEFORE THE EXECUTION OF THE TASK
 334   //
 335
 336   //slot #1
 337 //   OpenFile(1);
 338    fHmpHistList = new TList();
 339    fHmpHistList->SetOwner();
 340
 341    fHmpPesdPhmp = new TH2F("fHmpPesdPhmp","HMPID: ESD p - running p;HMP p (GeV/c);ESD p (Gev/c)",100,0,10,100,0,10);
 342    fHmpHistList->Add(fHmpPesdPhmp);
 343
 344    fHmpCkovPesd = new TH2F("fHmpCkovPesd","HMPID: ThetaCherenkov vs P;p_esd (GeV/c);#Theta_C;Entries",100,0,10,110,0,1.1);
 345    fHmpHistList->Add(fHmpCkovPesd);
 346
 347    fHmpCkovPhmp = new TH2F("fHmpCkovPhmp","HMPID: ThetaCherenkov vs P;p_hmp (GeV/c);#Theta_C;Entries",100,0,10,110,0,1.1);
 348    fHmpHistList->Add(fHmpCkovPhmp);
 349
 350    fHmpMipTrkDistX = new TH1F("fHmpMipTrkDistX","HMPID MIP-Track distance in local X;distance (cm);Entries",800,-20,20);
 351    fHmpHistList->Add(fHmpMipTrkDistX);
 352
 353    for (Int_t i=0; i<7; i++){
 354      TString title=Form("MIP-Track distance in local Y, Chamber %d;distance (cm);Entries",i);
 355      fHmpMipTrkDistPosY[i] = new TH1F(Form("fHmpMipTrkDistPosY%d",i),title.Data(),800,-20,20);
 356      fHmpHistList->Add(fHmpMipTrkDistPosY[i]);
 357      fHmpMipTrkDistNegY[i] = new TH1F(Form("fHmpMipTrkDistNegY%d",i),title.Data(),800,-20,20);
 358      fHmpHistList->Add(fHmpMipTrkDistNegY[i]);
 359    }
 360
 361    fHmpMipCharge3cm = new TH1F("fHmpMipCharge3cm","HMPID MIP Charge;MIP Charge (ADC);Entries",5001,-0.5,5000.5);
 362    fHmpHistList->Add(fHmpMipCharge3cm);
 363
 364    fHmpTrkFlags = new TH1F("fHmpTrkFlags","HMPID track flags",6,0,6);
 365    TString summary[6] =  {"NotPerformed","MipDistCut", "MipQdcCut", "NoPhotAccept", "kNoRad", "other"};
 366    for(Int_t ibin = 0; ibin < 6; ibin++) fHmpTrkFlags->GetXaxis()->SetBinLabel(ibin+1,Form("%i  %s",ibin+1,summary[ibin].Data()));
 367    fHmpHistList->Add(fHmpTrkFlags);
 368
 369    fPionEff = new TH1F("PionEff","Identified pions",fN1,0,fN1);
 370    fKaonEff = new TH1F("KaonEff","Identified kaons",fN1,0,fN1);
 371    fProtEff = new TH1F("ProtEff","Identified protons",fN2,0,fN2);
 372    fPionTot = new TH1I("PionTot","Total MC pions",fN1,0,fN1);
 373    fKaonTot = new TH1I("KaonTot","Total MC kaons",fN1,0,fN1);
 374    fProtTot = new TH1I("ProtTot","Total MC protons",fN2,0,fN2);
 375    fPionNot = new TH1F("PionNot","Misidentified pions",fN1,0,fN1);
 376    fKaonNot = new TH1F("KaonNot","Misidentified kaons",fN1,0,fN1);
 377    fProtNot = new TH1F("ProtNot","Misidentified protons",fN2,0,fN2);
 378    fPionCon = new TH1I("PionCon","Total not MC pions",fN1,0,fN1);
 379    fKaonCon = new TH1I("KaonCon","Total not MC kaons",fN1,0,fN1);
 380    fProtCon = new TH1I("ProtCon","Total not MC protons",fN2,0,fN2);
 381
 382    fHmpHistList->Add(fPionEff); fHmpHistList->Add(fKaonEff); fHmpHistList->Add(fProtEff);
 383    fHmpHistList->Add(fPionTot); fHmpHistList->Add(fKaonTot); fHmpHistList->Add(fProtTot);
 384    fHmpHistList->Add(fPionNot); fHmpHistList->Add(fKaonNot); fHmpHistList->Add(fProtNot);
 385    fHmpHistList->Add(fPionCon); fHmpHistList->Add(fKaonCon); fHmpHistList->Add(fProtCon);
 386
 387    PostData(1,fHmpHistList);
 388 }
 389
 390 //____________________________________________________________________________________________________________________________________
 391 Bool_t AliHMPIDTaskQA::Equal(Double_t x, Double_t y, Double_t tolerance)
 392 {
 393  return abs(x - y) <= tolerance ;
 394 }
 395
 396 #endif