ANALYSIS/AliAnalysisTaskPIDqa.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 /* $Id: AliAnalysisTaskPIDqa.cxx 43811 2010-09-23 14:13:31Z wiechula $ */
  17 #include <TList.h>
  18 #include <TVectorD.h>
  19 #include <TObjArray.h>
  20 #include <TH2.h>
  21 #include <TFile.h>
  22 #include <TPRegexp.h>
  23 #include <TChain.h>
  24 #include <TF1.h>
  25 #include <TSpline.h>
  26
  27 #include <AliAnalysisManager.h>
  28 #include <AliInputEventHandler.h>
  29 #include <AliVEventHandler.h>
  30 #include <AliVEvent.h>
  31 #include <AliVParticle.h>
  32 #include <AliVTrack.h>
  33 #include <AliLog.h>
  34 #include <AliPID.h>
  35 #include <AliPIDResponse.h>
  36 #include <AliITSPIDResponse.h>
  37 #include <AliTPCPIDResponse.h>
  38 #include <AliTRDPIDResponse.h>
  39 #include <AliTOFPIDResponse.h>
  40 #include <AliTPCdEdxInfo.h>
  41
  42 #include <AliESDEvent.h>
  43 #include <AliAODEvent.h>
  44 #include <AliESDv0.h>
  45 #include <AliAODv0.h>
  46 #include <AliESDv0KineCuts.h>
  47 #include <AliESDtrackCuts.h>
  48
  49 #include <AliMCEvent.h>
  50
  51 #include "AliAnalysisTaskPIDqa.h"
  52
  53
  54 ClassImp(AliAnalysisTaskPIDqa)
  55
  56 //______________________________________________________________________________
  57 AliAnalysisTaskPIDqa::AliAnalysisTaskPIDqa():
  58 AliAnalysisTaskSE(),
  59 fPIDResponse(0x0),
  60 fV0cuts(0x0),
  61 fV0electrons(0x0),
  62 fV0pions(0x0),
  63 fV0kaons(0x0),
  64 fV0protons(0x0),
  65 fListQA(0x0),
  66 fListQAits(0x0),
  67 fListQAitsSA(0x0),
  68 fListQAitsPureSA(0x0),
  69 fListQAtpc(0x0),
  70 fListQAtpcBasic(0x0),
  71 fListQAtpcMCtruth(0x0),
  72 fListQAtpcHybrid(0x0),
  73 fListQAtpcOROChigh(0x0),
  74 fListQAtpcV0(0x0),
  75 fListQAtrd(0x0),
  76 fListQAtrdNsig(0x0),
  77 fListQAtrdNsigTPCTOF(0x0),
  78 fListQAtof(0x0),
  79 fListQAt0(0x0),
  80 fListQAemcal(0x0),
  81 fListQAhmpid(0x0),
  82 fListQAtofhmpid(0x0),
  83 fListQAtpctof(0x0),
  84 fListQAV0(0x0),
  85 fListQAinfo(0x0)
  86 {
  87   //
  88   // Dummy constructor
  89   //
  90 }
  91
  92 //______________________________________________________________________________
  93 AliAnalysisTaskPIDqa::AliAnalysisTaskPIDqa(const char* name):
  94 AliAnalysisTaskSE(name),
  95 fPIDResponse(0x0),
  96 fV0cuts(0x0),
  97 fV0electrons(0x0),
  98 fV0pions(0x0),
  99 fV0kaons(0x0),
 100 fV0protons(0x0),
 101 fListQA(0x0),
 102 fListQAits(0x0),
 103 fListQAitsSA(0x0),
 104 fListQAitsPureSA(0x0),
 105 fListQAtpc(0x0),
 106 fListQAtpcBasic(0x0),
 107 fListQAtpcMCtruth(0x0),
 108 fListQAtpcHybrid(0x0),
 109 fListQAtpcOROChigh(0x0),
 110 fListQAtpcV0(0x0),
 111 fListQAtrd(0x0),
 112 fListQAtrdNsig(0x0),
 113 fListQAtrdNsigTPCTOF(0x0),
 114 fListQAtof(0x0),
 115 fListQAt0(0x0),
 116 fListQAemcal(0x0),
 117 fListQAhmpid(0x0),
 118 fListQAtofhmpid(0x0),
 119 fListQAtpctof(0x0),
 120 fListQAV0(0x0),
 121 fListQAinfo(0x0)
 122 {
 123   //
 124   // Default constructor
 125   //
 126   DefineInput(0,TChain::Class());
 127   DefineOutput(1,TList::Class());
 128 }
 129
 130 //______________________________________________________________________________
 131 AliAnalysisTaskPIDqa::~AliAnalysisTaskPIDqa()
 132 {
 133   //
 134   // Destructor
 135   //
 136
 137   delete fV0cuts;
 138   delete fV0electrons;
 139   delete fV0pions;
 140   delete fV0kaons;
 141   delete fV0protons;
 142
 143   if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fListQA;
 144 }
 145
 146 //______________________________________________________________________________
 147 void AliAnalysisTaskPIDqa::UserCreateOutputObjects()
 148 {
 149   //
 150   // Create the output QA objects
 151   //
 152
 153   AliLog::SetClassDebugLevel("AliAnalysisTaskPIDqa",10);
 154
 155   //input hander
 156   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 157   AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
 158   if (!inputHandler) AliFatal("Input handler needed");
 159
 160   //pid response object
 161   fPIDResponse=inputHandler->GetPIDResponse();
 162   if (!fPIDResponse) AliError("PIDResponse object was not created");
 163
 164   // V0 Kine cuts
 165   fV0cuts = new AliESDv0KineCuts;
 166
 167   // V0 PID Obj arrays
 168   fV0electrons = new TObjArray;
 169   fV0pions     = new TObjArray;
 170   fV0kaons     = new TObjArray;
 171   fV0protons   = new TObjArray;
 172
 173   //
 174   fListQA=new TList;
 175   fListQA->SetOwner();
 176
 177   fListQAits=new TList;
 178   fListQAits->SetOwner();
 179   fListQAits->SetName("ITS");
 180
 181   fListQAitsSA=new TList;
 182   fListQAitsSA->SetOwner();
 183   fListQAitsSA->SetName("ITS_SA");
 184
 185   fListQAitsPureSA=new TList;
 186   fListQAitsPureSA->SetOwner();
 187   fListQAitsPureSA->SetName("ITS_PureSA");
 188
 189   fListQAtpc=new TList;
 190   fListQAtpc->SetOwner();
 191   fListQAtpc->SetName("TPC");
 192
 193   fListQAtrd=new TList;
 194   fListQAtrd->SetOwner();
 195   fListQAtrd->SetName("TRD");
 196
 197   fListQAtrdNsig=new TList;
 198   fListQAtrdNsig->SetOwner();
 199   fListQAtrdNsig->SetName("TRDnSigma");
 200
 201   fListQAtrdNsigTPCTOF=new TList;
 202   fListQAtrdNsigTPCTOF->SetOwner();
 203   fListQAtrdNsigTPCTOF->SetName("TRDnSigma_TPCTOF");
 204
 205   fListQAtof=new TList;
 206   fListQAtof->SetOwner();
 207   fListQAtof->SetName("TOF");
 208
 209   fListQAt0=new TList;
 210   fListQAt0->SetOwner();
 211   fListQAt0->SetName("T0");
 212
 213   fListQAemcal=new TList;
 214   fListQAemcal->SetOwner();
 215   fListQAemcal->SetName("EMCAL");
 216
 217   fListQAhmpid=new TList;
 218   fListQAhmpid->SetOwner();
 219   fListQAhmpid->SetName("HMPID");
 220
 221   fListQAtpctof=new TList;
 222   fListQAtpctof->SetOwner();
 223   fListQAtpctof->SetName("TPC_TOF");
 224
 225   fListQAtofhmpid=new TList;
 226   fListQAtofhmpid->SetOwner();
 227   fListQAtofhmpid->SetName("TOF_HMPID");
 228
 229   fListQAV0=new TList;
 230   fListQAV0->SetOwner();
 231   fListQAV0->SetName("V0decay");
 232
 233   fListQAinfo=new TList;
 234   fListQAinfo->SetOwner();
 235   fListQAinfo->SetName("QAinfo");
 236
 237   fListQA->Add(fListQAits);
 238   fListQA->Add(fListQAitsSA);
 239   fListQA->Add(fListQAitsPureSA);
 240   fListQA->Add(fListQAtpc);
 241   fListQA->Add(fListQAtrd);
 242   fListQA->Add(fListQAtof);
 243   fListQA->Add(fListQAt0);
 244   fListQA->Add(fListQAemcal);
 245   fListQA->Add(fListQAhmpid);
 246   fListQA->Add(fListQAtpctof);
 247   fListQA->Add(fListQAtofhmpid);
 248   fListQA->Add(fListQAV0);
 249   fListQA->Add(fListQAinfo);
 250
 251   SetupITSqa();
 252 //  SetupTPCqa(kFALSE, kTRUE, kFALSE);
 253   SetupTRDqa();
 254   SetupTOFqa();
 255   SetupT0qa();
 256   SetupEMCALqa();
 257   SetupHMPIDqa();
 258   SetupTPCTOFqa();
 259   SetupTOFHMPIDqa();
 260   SetupV0qa();
 261   SetupQAinfo();
 262
 263   PostData(1,fListQA);
 264 }
 265
 266
 267 //______________________________________________________________________________
 268 void AliAnalysisTaskPIDqa::UserExec(Option_t */*option*/)
 269 {
 270   //
 271   // Setup the PID response functions and fill the QA histograms
 272   //
 273
 274   AliVEvent *event=InputEvent();
 275   if (!event||!fPIDResponse) return;
 276
 277   // Start with the V0 task (only possible for ESDs?)
 278   FillV0PIDlist();
 279
 280   FillITSqa();
 281   FillTPCqa();
 282   FillTRDqa();
 283   FillTOFqa();
 284   FillEMCALqa();
 285   FillHMPIDqa();
 286   FillT0qa();
 287
 288   //combined detector QA
 289   FillTPCTOFqa();
 290   FillTOFHMPIDqa();
 291
 292   // Clear the V0 PID arrays
 293   ClearV0PIDlist();
 294
 295   //QA info
 296   FillQAinfo();
 297
 298   PostData(1,fListQA);
 299 }
 300
 301 //______________________________________________________________________________
 302 void  AliAnalysisTaskPIDqa::FillV0PIDlist(){
 303
 304   //
 305   // Fill the PID object arrays holding the pointers to identified particle tracks
 306   //
 307
 308   // Dynamic cast to ESD events (DO NOTHING for AOD events)
 309   AliESDEvent *event = dynamic_cast<AliESDEvent *>(InputEvent());
 310   if ( !event )  return;
 311
 312   if(TString(event->GetBeamType())=="Pb-Pb" || TString(event->GetBeamType())=="A-A"){
 313     fV0cuts->SetMode(AliESDv0KineCuts::kPurity,AliESDv0KineCuts::kPbPb);
 314   }
 315   else{
 316     fV0cuts->SetMode(AliESDv0KineCuts::kPurity,AliESDv0KineCuts::kPP);
 317   }
 318
 319   // V0 selection
 320   // set event
 321   fV0cuts->SetEvent(event);
 322
 323   // loop over V0 particles
 324   for(Int_t iv0=0; iv0<event->GetNumberOfV0s();iv0++){
 325
 326     AliESDv0 *v0 = (AliESDv0 *) event->GetV0(iv0);
 327
 328     if(!v0) continue;
 329     if(v0->GetOnFlyStatus()) continue;
 330
 331     // Get the particle selection
 332     Bool_t foundV0 = kFALSE;
 333     Int_t pdgV0, pdgP, pdgN;
 334
 335     foundV0 = fV0cuts->ProcessV0(v0, pdgV0, pdgP, pdgN);
 336     if(!foundV0) continue;
 337
 338     Int_t iTrackP = v0->GetPindex();  // positive track
 339     Int_t iTrackN = v0->GetNindex();  // negative track
 340
 341     // v0 Armenteros plot (QA)
 342     Float_t armVar[2] = {0.0,0.0};
 343     fV0cuts->Armenteros(v0, armVar);
 344
 345     TH2 *h=(TH2*)fListQAV0->At(0);
 346     if (!h) continue;
 347     h->Fill(armVar[0],armVar[1]);
 348
 349     // fill the Object arrays
 350     // positive particles
 351     if( pdgP == -11){
 352       fV0electrons->Add((AliVTrack*)event->GetTrack(iTrackP));
 353     }
 354     else if( pdgP == 211){
 355       fV0pions->Add((AliVTrack*)event->GetTrack(iTrackP));
 356     }
 357     else if( pdgP == 321){
 358       fV0kaons->Add((AliVTrack*)event->GetTrack(iTrackP));
 359     }
 360     else if( pdgP == 2212){
 361       fV0protons->Add((AliVTrack*)event->GetTrack(iTrackP));
 362     }
 363
 364     // negative particles
 365     if( pdgN == 11){
 366       fV0electrons->Add((AliVTrack*)event->GetTrack(iTrackN));
 367     }
 368     else if( pdgN == -211){
 369       fV0pions->Add((AliVTrack*)event->GetTrack(iTrackN));
 370     }
 371     else if( pdgN == -321){
 372       fV0kaons->Add((AliVTrack*)event->GetTrack(iTrackN));
 373     }
 374     else if( pdgN == -2212){
 375       fV0protons->Add((AliVTrack*)event->GetTrack(iTrackN));
 376     }
 377
 378
 379   }
 380 }
 381 //______________________________________________________________________________
 382 void  AliAnalysisTaskPIDqa::ClearV0PIDlist(){
 383
 384   //
 385   // Clear the PID object arrays
 386   //
 387
 388   fV0electrons->Clear();
 389   fV0pions->Clear();
 390   fV0kaons->Clear();
 391   fV0protons->Clear();
 392
 393 }
 394 //______________________________________________________________________________
 395 void AliAnalysisTaskPIDqa::FillITSqa()
 396 {
 397   //
 398   // Fill PID qa histograms for the ITS
 399   //
 400
 401   AliVEvent *event=InputEvent();
 402
 403   Int_t ntracks=event->GetNumberOfTracks();
 404   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 405     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 406     ULong_t status=track->GetStatus();
 407     // not that nice. status bits not in virtual interface
 408     // ITS refit + ITS pid selection
 409     if (!( ( (status & AliVTrack::kITSrefit)==AliVTrack::kITSrefit ) ||
 410            ! ( (status & AliVTrack::kITSpid  )==AliVTrack::kITSpid   ) )) continue;
 411     Double_t mom=track->P();
 412
 413     TList *theList = 0x0;
 414     if(( (status & AliVTrack::kTPCin)==AliVTrack::kTPCin )){
 415       //ITS+TPC tracks
 416       theList=fListQAits;
 417     }else{
 418       if(!( (status & AliVTrack::kITSpureSA)==AliVTrack::kITSpureSA )){
 419         //ITS Standalone tracks
 420         theList=fListQAitsSA;
 421       }else{
 422         //ITS Pure Standalone tracks
 423         theList=fListQAitsPureSA;
 424       }
 425     }
 426
 427
 428     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
 429       TH2 *h=(TH2*)theList->At(ispecie);
 430       if (!h) continue;
 431       Double_t nSigma=fPIDResponse->NumberOfSigmasITS(track, (AliPID::EParticleType)ispecie);
 432       h->Fill(mom,nSigma);
 433     }
 434     TH2 *h=(TH2*)theList->At(AliPID::kSPECIESC);
 435     if (h) {
 436       Double_t sig=track->GetITSsignal();
 437       h->Fill(mom,sig);
 438     }
 439   }
 440 }
 441
 442
 443 //______________________________________________________________________________
 444 void AliAnalysisTaskPIDqa::FillTPCHistogramsSignal(TList *sublist, Int_t scenario, AliVTrack *track, Int_t mult)
 445 {
 446   //
 447   // Fill PID qa histograms for the TPC: Fill the histograms for the TPC signal for different settings
 448   //
 449
 450   AliMCEvent *eventMC=MCEvent();  // MC event for MC truth PID
 451
 452   Double_t mom=0.;      // track momentum
 453   Double_t eta=0.;      // track eta
 454   Double_t sig=0.;      // TPC dE/dx signal
 455   Double_t sigStd=0.;         // TPC dE/dx signal (standard = all ROCs)
 456   Double_t sigIROC=0.;        // TPC dE/dx signal (IROC)
 457   Double_t sigOROCmedium=0.;  // TPC dE/dx signal (OROCmedium)
 458   Double_t sigOROClong=0.;    // TPC dE/dx signal (OROClong)
 459   Double_t eleLineDist=0.;    // difference between TPC signal and electron expectation
 460   Int_t trackLabel=0;   // label of the AliVTrack to identify the corresponding MCtrack
 461   Int_t pdgCode=0;      // pdgcode of MC track for MC truth scenario
 462   Int_t pdgCodeAbs=0;   // absolute value of pdgcode to get both particles and antiparticles
 463   Int_t iSigMax=1;      // number of TPC signals (std = 1, set automatically higher if available)
 464   Int_t nSpecies=0;     // number of particle species under study
 465   Int_t count=0;        // counter for the number of plot sets for all species (i.e. nsigma vs. p, eta and mult)
 466
 467   mom=track->GetTPCmomentum();
 468   eta=track->Eta();
 469   sigStd=track->GetTPCsignal();
 470
 471   eleLineDist=sigStd-fPIDResponse->GetTPCResponse().GetExpectedSignal(track,AliPID::kElectron);
 472
 473   // Get number of particle species (less for V0 candidates = scenarios 40-44)
 474   if (scenario > 39) nSpecies=(Int_t)AliPID::kSPECIES;
 475   else nSpecies=(Int_t)AliPID::kSPECIESC;
 476
 477   // Set number of plot sets for all species
 478   // (i.e. only nsigma vs. p => count=1; also vs. eta and mult => count=3)
 479   if ( scenario == 1 || scenario > 39) count=3;
 480   else count=1;
 481
 482   // Get MC track ( --> can be deleted if TPC signal is NOT filled for scenario=1 (MC truth)
 483   if (eventMC) {
 484     trackLabel=TMath::Abs(track->GetLabel());
 485     AliVTrack *mcTrack=(AliVTrack*)eventMC->GetTrack(trackLabel);
 486     pdgCode=mcTrack->PdgCode();
 487     pdgCodeAbs=TMath::Abs(pdgCode);
 488   }
 489
 490   // Get TPC dE/dx info and different TPC signals (IROC, OROCmedium, OROClong)
 491   AliTPCdEdxInfo* fTPCdEdxInfo = 0x0;
 492   fTPCdEdxInfo = track->GetTPCdEdxInfo();
 493
 494   if (fTPCdEdxInfo) {
 495     sigIROC=fTPCdEdxInfo->GetTPCsignalShortPad();
 496     sigOROCmedium=fTPCdEdxInfo->GetTPCsignalMediumPad();
 497     sigOROClong=fTPCdEdxInfo->GetTPCsignalLongPad();
 498     iSigMax=4;
 499
 500     //printf("mom = %.3f  sigStd = %.3f  sigIROC = %.3f  sigOROCmedium = %.3f  sigOROClong = %.3f \n",mom,sigStd,sigIROC,sigOROCmedium,sigOROClong);
 501   }
 502
 503
 504   // TPC signal for all particles vs. momentum (standard, IROC, OROCmedium, OROClong)
 505   TH2 *h1std=(TH2*)sublist->At(count*nSpecies+4);
 506   if (h1std) {
 507     h1std->Fill(mom,sigStd);
 508   }
 509
 510   TH2 *h1iroc=(TH2*)sublist->At(count*nSpecies+5);
 511   if ( h1iroc && sigIROC ) {
 512     h1iroc->Fill(mom,sigIROC);
 513   }
 514
 515   TH2 *h1orocm=(TH2*)sublist->At(count*nSpecies+6);
 516   if  (h1orocm && sigOROCmedium ) {
 517     h1orocm->Fill(mom,sigOROCmedium);
 518   }
 519
 520   TH2 *h1orocl=(TH2*)sublist->At(count*nSpecies+7);
 521   if ( h1orocl && sigOROClong ) {
 522     h1orocl->Fill(mom,sigOROClong);
 523   }
 524
 525
 526   // - Beginn: MIP pions: TPC signal vs. eta, TPC signal vs. mult -
 527   if (mom>0.45 && mom<0.5 && sigStd>40 && sigStd<60) {
 528
 529     Bool_t isPionMC=kTRUE;
 530
 531     if (scenario == 1) {
 532       if ( pdgCodeAbs != 211 && pdgCodeAbs != 111 ) isPionMC=kFALSE;
 533     }
 534
 535     // MIP pions: TPC signal vs. eta (standard, IROC, OROCmedium, OROClong)
 536     for (Int_t iSig=0; iSig<iSigMax; iSig++) {
 537       if (iSig==0) sig=sigStd;
 538       else if (iSig==1) sig=sigIROC;
 539       else if (iSig==2) sig=sigOROCmedium;
 540       else if (iSig==3) sig=sigOROClong;
 541
 542       TH2 *h2=(TH2*)sublist->At(count*nSpecies+8+iSig);
 543       if ( h2 && isPionMC ) {
 544         h2->Fill(eta,sig);
 545       }
 546     }
 547
 548     // MIP pions: TPC signal vs. mult (standard, IROC, OROCmedium, OROClong)
 549     for (Int_t iSig=0; iSig<iSigMax; iSig++) {
 550       if (iSig==0) sig=sigStd;
 551       else if (iSig==1) sig=sigIROC;
 552       else if (iSig==2) sig=sigOROCmedium;
 553       else if (iSig==3) sig=sigOROClong;
 554
 555       TH2 *h3=(TH2*)sublist->At(count*nSpecies+12+iSig);
 556       if ( h3 && isPionMC && mult > 0 ) {
 557         h3->Fill(mult,sig);
 558       }
 559     }
 560   } // - End: MIP pions -
 561
 562   // - Beginn: Electrons: TPC signal vs. eta, TPC signal vs. mult -
 563   if (mom>0.32 && mom<0.38 && eleLineDist>-10. && eleLineDist<15.) {
 564
 565     Bool_t isElectronMC=kTRUE;
 566
 567     if (scenario == 1) {
 568       if ( pdgCodeAbs != 11 ) isElectronMC=kFALSE;
 569     }
 570
 571     // Electrons: TPC signal vs. eta (standard, IROC, OROCmedium, OROClong)
 572     for (Int_t iSig=0; iSig<iSigMax; iSig++) {
 573       if (iSig==0) sig=sigStd;
 574       else if (iSig==1) sig=sigIROC;
 575       else if (iSig==2) sig=sigOROCmedium;
 576       else if (iSig==3) sig=sigOROClong;
 577
 578       TH2 *h4=(TH2*)sublist->At(count*nSpecies+16+iSig);
 579       if ( h4 && isElectronMC ) {
 580         h4->Fill(eta,sig);
 581       }
 582     }
 583
 584     // Electrons: TPC signal vs. mult (standard, IROC, OROCmedium, OROClong)
 585     for (Int_t iSig=0; iSig<iSigMax; iSig++) {
 586       if (iSig==0) sig=sigStd;
 587       else if (iSig==1) sig=sigIROC;
 588       else if (iSig==2) sig=sigOROCmedium;
 589       else if (iSig==3) sig=sigOROClong;
 590
 591       TH2 *h5=(TH2*)sublist->At(count*nSpecies+20+iSig);
 592       if ( h5 && isElectronMC && mult > 0 ) {
 593         h5->Fill(mult,sig);
 594       }
 595     }
 596   } // - End: Electrons -
 597
 598 }
 599
 600 //______________________________________________________________________________
 601 void AliAnalysisTaskPIDqa::FillTPCHistogramsNsigma(TList *sublist, Int_t scenario, AliVTrack *track, Int_t mult)
 602 {
 603   //
 604   // Fill PID qa histograms for the TPC: Fill the histograms for TPC Nsigma for different settings
 605   //
 606
 607   AliMCEvent *eventMC=MCEvent();  // MC event for MC truth PID
 608
 609   Double_t mom=0.;      // track momentum
 610   Double_t eta=0.;      // track eta
 611   Double_t nSigma=0.;   // number of sigmas wrt. expected signal
 612   Double_t sig=0.;      // TPC dE/dx signal
 613   Double_t eleLineDist=0.;  // difference between TPC signal and electron expectation
 614   Int_t trackLabel=0;   // label of the AliVTrack to identify the corresponding MCtrack
 615   Int_t pdgCode=0;      // pdgcode of MC track for MC truth scenario
 616   Int_t pdgCodeAbs=0;   // absolute value of pdgcode to get both particles and antiparticles
 617   Int_t nSpecies=0;     // number of particle species under study
 618   Int_t count=0;        // counter for the number of plot sets for all species (i.e. vs. p, eta and mult)
 619
 620   mom=track->GetTPCmomentum();
 621   eta=track->Eta();
 622   sig=track->GetTPCsignal();
 623
 624   eleLineDist=sig-fPIDResponse->GetTPCResponse().GetExpectedSignal(track,AliPID::kElectron);
 625
 626   // Get number of particle species (less for V0 candidates = scenarios 40-44)
 627   if (scenario > 39) nSpecies=(Int_t)AliPID::kSPECIES;
 628   else nSpecies=(Int_t)AliPID::kSPECIESC;
 629
 630   // Set number of plot sets for all species
 631   // (i.e. only vs. p => count=1; also vs. eta and mult => count=3)
 632   if ( scenario == 1 || scenario > 39 ) count=3;
 633   else count=1;
 634
 635   // Get MC track
 636   if (eventMC) {
 637     trackLabel=TMath::Abs(track->GetLabel());
 638     AliVTrack *mcTrack=(AliVTrack*)eventMC->GetTrack(trackLabel);
 639     pdgCode=mcTrack->PdgCode();
 640     pdgCodeAbs=TMath::Abs(pdgCode);
 641   }
 642
 643
 644   // - Beginn: Nsigma vs. p, vs. eta and vs. multiplicity for different particle species -
 645   for (Int_t ispecie=0; ispecie<nSpecies; ++ispecie){
 646
 647     TH2 *h=(TH2*)sublist->At(ispecie);
 648     if (!h) continue;
 649
 650     if (scenario == 1) {
 651       if ( ispecie == 0 && pdgCodeAbs != 11 ) continue;  // Electron
 652       if ( ispecie == 1 && pdgCodeAbs != 13 ) continue;  // Muon
 653       if ( ispecie == 2 && pdgCodeAbs != 211 && pdgCodeAbs!=111 ) continue;  // Pion
 654       if ( ispecie == 3 && pdgCodeAbs != 321 && pdgCodeAbs!=311 ) continue;  // Kaon
 655       if ( ispecie == 4 && pdgCodeAbs != 2212 ) continue;  // Proton
 656       if ( ispecie == 5 && pdgCodeAbs != 1000010020 ) continue;  // Deuteron
 657       if ( ispecie == 6 && pdgCodeAbs != 1000010030 ) continue;  // Triton
 658       if ( ispecie == 7 && pdgCodeAbs != 1000020030 ) continue;  // Helium-3
 659       if ( ispecie == 8 && pdgCodeAbs != 1000020040 ) continue;  // Alpha
 660     }
 661     else if (scenario > 39) {
 662       if ( ispecie == 0 && scenario != 40 ) continue;  // Electron
 663       if ( ispecie == 1 ) continue;  // Muon
 664       if ( ispecie == 2 && scenario != 42 ) continue;  // Pion
 665       if ( ispecie == 3 && scenario != 43 ) continue;  // Kaon
 666       if ( ispecie == 4 && scenario != 44 ) continue;  // Proton
 667     }
 668
 669     if (scenario == 2) {
 670       nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie, AliTPCPIDResponse::kdEdxHybrid);
 671     }
 672     else if (scenario == 3) {
 673       nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie, AliTPCPIDResponse::kdEdxOROC);
 674     }
 675     else {
 676       nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
 677     }
 678
 679     h->Fill(mom,nSigma);
 680
 681     if (count == 3) {
 682       TH2 *hEta=(TH2*)sublist->At(ispecie+nSpecies);
 683       TH2 *hMult=(TH2*)sublist->At(ispecie+2*nSpecies);
 684
 685       if ( hEta ) hEta->Fill(eta,nSigma);
 686       if ( hMult && mult > 0 ) hMult->Fill(mult,nSigma);
 687     }
 688   } // - End: different particle species -
 689
 690
 691   // -- Beginn: Fill histograms for MIP pions and electrons (only for some scenarios) --
 692   if ( scenario == 0 || scenario == 2 || scenario == 3 ) {
 693
 694     // - Beginn: MIP pions: Nsigma vs. eta, Nsigma vs. mult -
 695     if (mom>0.45 && mom<0.5 && sig>40 && sig<60) {
 696
 697       Bool_t isPionMC=kTRUE;
 698
 699       TH2 *h1=(TH2*)sublist->At(count*nSpecies);
 700       if (h1) {
 701         if (scenario == 1) {
 702           if ( pdgCodeAbs != 211 && pdgCodeAbs != 111 ) isPionMC=kFALSE;
 703           if (isPionMC) {
 704             nSigma=fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion);
 705           }
 706         }
 707         else if (scenario == 2) {
 708           nSigma=fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion, AliTPCPIDResponse::kdEdxHybrid);
 709         }
 710         else if (scenario == 3) {
 711           nSigma=fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion, AliTPCPIDResponse::kdEdxOROC);
 712         }
 713         else nSigma=fPIDResponse->NumberOfSigmasTPC(track, AliPID::kPion);
 714
 715         if (isPionMC) h1->Fill(eta,nSigma);
 716       }
 717
 718       TH2 *h2m=(TH2*)sublist->At(count*nSpecies+1);
 719       if ( h2m && isPionMC && mult > 0 ) {
 720         h2m->Fill(mult,nSigma);
 721       }
 722
 723     } // - End: MIP pions -
 724
 725     // - Beginn: Electrons: Nsigma vs. eta, Nsigma vs. mult -
 726     if (mom>0.32 && mom<0.38 && eleLineDist>-10. && eleLineDist<15.) {
 727
 728       Bool_t isElectronMC=kTRUE;
 729
 730       TH2 *h3=(TH2*)sublist->At(count*nSpecies+2);
 731       if (h3) {
 732         if (scenario == 1) {
 733           if ( pdgCodeAbs != 11 ) isElectronMC=kFALSE;
 734           if (isElectronMC) {
 735             nSigma=fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
 736           }
 737         }
 738         if (scenario == 2) {
 739           nSigma=fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxHybrid);
 740         }
 741         else if (scenario == 3) {
 742           nSigma=fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron, AliTPCPIDResponse::kdEdxOROC);
 743         }
 744         else nSigma=fPIDResponse->NumberOfSigmasTPC(track, AliPID::kElectron);
 745
 746         if (isElectronMC) h3->Fill(eta,nSigma);
 747       }
 748
 749       TH2 *h4m=(TH2*)sublist->At(count*nSpecies+3);
 750       if ( h4m && isElectronMC && mult > 0 ) {
 751         h4m->Fill(mult,nSigma);
 752       }
 753
 754     } // - End: Electrons -
 755   } // -- End: Fill histograms for MIP pions and electrons --
 756
 757 }
 758
 759 //______________________________________________________________________________
 760 void AliAnalysisTaskPIDqa::FillTPCqa()
 761 {
 762   //
 763   // Fill PID qa histograms for the TPC
 764   //
 765
 766   // switches for the different scenarios
 767   Bool_t scBasic=1;     // default/basic
 768   Bool_t scMCtruth=1;   // for MC truth tracks
 769   Bool_t scHybrid=1;    // for hybrid PID (only LHC11h)
 770   Bool_t scOROChigh=1;  // only OROC signal (only LHC11h)
 771   Bool_t scV0=1;        // for V0 candidates (only for ESDs available)
 772   Int_t scCounter=0;    // counter of scenarios, used for the histograms at the end of FillTPCqa
 773
 774   // input handler
 775   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 776   AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
 777   if (!inputHandler) AliFatal("Input handler needed");
 778
 779   AliVEvent *event=InputEvent();
 780
 781   // ESD or AOD event needed to get reference multiplicity (not in AliVEvent)
 782   AliAODEvent *fAODevent = 0x0;   // AOD event
 783   AliESDEvent *fESDevent = 0x0;   // ESD event
 784   AliESDtrackCuts *esdTrackCuts = 0x0;  // ESD track Cuts (ref mult is in AliESDtrackCuts)
 785
 786   Double_t eta=0.;    // track eta
 787   Int_t mult=0;       // event multiplicity (TPConlyRefMult)
 788   //Int_t nacc=0;       // counter for accepted multiplicity
 789
 790   // Check for MC
 791   scMCtruth=(MCEvent()!=0x0);
 792
 793   // Check if period is data LHC11h by checking if
 794   // the splines for ALLhigh have been set by AliPIDResponse
 795   AliTPCPIDResponse &tpcResp=fPIDResponse->GetTPCResponse();
 796   if (tpcResp.GetResponseFunction(AliPID::kPion, AliTPCPIDResponse::kALLhigh)==0x0) {
 797     scHybrid   = kFALSE;
 798     scOROChigh = kFALSE;
 799   }
 800
 801   // Check if "ESD" or "AOD" and get the corresponding event and the beam type (or centrality)
 802   TString analysisType = inputHandler->GetDataType(); // can be "ESD" or "AOD"
 803   if (analysisType == "ESD") {
 804     fESDevent = dynamic_cast<AliESDEvent*>( InputEvent() );
 805     esdTrackCuts = new AliESDtrackCuts("esdTrackCuts");
 806     //printf("\n--- New event - event type = ESD \n");
 807   }
 808   else if (analysisType == "AOD") {
 809     fAODevent = dynamic_cast<AliAODEvent*>( InputEvent() );
 810     //printf("\n--- New event - event type = AOD \n");
 811
 812     // disable V0 scenario, because V0s are not available for AODs in this current implementation
 813     scV0=0;
 814   }
 815
 816   // Check if Basic list is already created
 817   // If not: Go to SetupTPCqa and creat lists and histograms
 818   if(!fListQAtpcBasic) {
 819     //printf("\n--- No list QA TPC Basic found -> go to SetupTPCqa! ---\n");
 820     SetupTPCqa(scMCtruth, scHybrid, scV0);
 821   }
 822
 823   // Get the number of scenarios by counting those, which are switched on
 824   if (scBasic) scCounter++;
 825   if (scMCtruth) scCounter++;
 826   if (scHybrid) scCounter++;
 827   if (scOROChigh) scCounter++;
 828   if (scV0) scCounter++;
 829
 830   // Get reference multiplicity for ESDs
 831   if ( analysisType == "ESD" && esdTrackCuts ) {
 832     mult=esdTrackCuts->GetReferenceMultiplicity(fESDevent,kTRUE);
 833   }
 834
 835   // Get reference multiplicity for AODs
 836   if ( analysisType == "AOD" && fAODevent ) {
 837     mult=fAODevent->GetHeader()->GetTPConlyRefMultiplicity();
 838   }
 839
 840   /*if (mult < 0) {
 841     printf("Reference multiplicity not available \n");
 842     //return;
 843   }*/
 844
 845   //printf("The multiplicity is = %i ",mult);
 846
 847
 848   // -- Begin: track loop --
 849   Int_t ntracks=event->GetNumberOfTracks();
 850   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 851     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 852
 853     //
 854     //basic track cuts
 855     //
 856     ULong_t status=track->GetStatus();
 857     // not that nice. status bits not in virtual interface
 858     // TPC refit + ITS refit + TPC pid
 859     if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 860         !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
 861
 862     // The TPC pid cut removes the light nuclei (>5 sigma from proton line)
 863     //||        !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  )
 864     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 865     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 866     if (track->GetTPCNclsF()>0) {
 867       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 868     }
 869
 870     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 871
 872     eta=track->Eta();
 873     if ( TMath::Abs(eta)>0.9 ) continue;
 874
 875     //nacc++; // counter for accepted multiplicity
 876
 877     // the default ("basic") scenario
 878     if (scBasic == 1) {
 879       FillTPCHistogramsNsigma(fListQAtpcBasic,0,track,mult);
 880       FillTPCHistogramsSignal(fListQAtpcBasic,0,track,mult);
 881     }
 882
 883     // only MC truth identified particles
 884     if (scMCtruth == 1) {
 885       FillTPCHistogramsNsigma(fListQAtpcMCtruth,1,track,mult);
 886     }
 887
 888     // the "hybrid" scenario (only for LHC11h)
 889     if (scHybrid == 1) {
 890       FillTPCHistogramsNsigma(fListQAtpcHybrid,2,track,mult);
 891     }
 892
 893     // the "OROC high" scenario (only for LHC11h)
 894     if (scOROChigh == 1) {
 895       FillTPCHistogramsNsigma(fListQAtpcOROChigh,3,track,mult);
 896     }
 897
 898   } // -- End: track loop --
 899
 900
 901   // -- Begin: track loops for V0 candidates --
 902   if (scV0 == 1) {
 903
 904     // - Begin: track loop for electrons from V0 -
 905     for(Int_t itrack = 0; itrack < fV0electrons->GetEntries(); itrack++){
 906       AliVTrack *track=(AliVTrack*)fV0electrons->At(itrack);
 907
 908       //
 909       //basic track cuts
 910       //
 911       ULong_t status=track->GetStatus();
 912       // not that nice. status bits not in virtual interface
 913       // TPC refit + ITS refit + TPC pid
 914       if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 915           !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
 916
 917       // The TPC pid cut removes the light nuclei (>5 sigma from proton line)
 918       //||        !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  )
 919       Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 920       Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 921       if (track->GetTPCNclsF()>0) {
 922         ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 923       }
 924
 925       if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 926
 927       eta=track->Eta();
 928       if ( TMath::Abs(eta)>0.9 ) continue;
 929
 930       // fill histograms for V0 candidates
 931       FillTPCHistogramsNsigma(fListQAtpcV0,40,track,mult);
 932
 933     } // - End: track loop for electrons from V0 -
 934
 935
 936     // - Begin: track loop for pions from V0 -
 937     for(Int_t itrack = 0; itrack < fV0pions->GetEntries(); itrack++){
 938       AliVTrack *track=(AliVTrack*)fV0pions->At(itrack);
 939
 940       //
 941       //basic track cuts
 942       //
 943       ULong_t status=track->GetStatus();
 944       // not that nice. status bits not in virtual interface
 945       // TPC refit + ITS refit + TPC pid
 946       if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 947           !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
 948
 949       // The TPC pid cut removes the light nuclei (>5 sigma from proton line)
 950       //||        !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  )
 951       Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 952       Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 953       if (track->GetTPCNclsF()>0) {
 954         ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 955       }
 956
 957       if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 958
 959       eta=track->Eta();
 960       if ( TMath::Abs(eta)>0.9 ) continue;
 961
 962       // fill histograms for V0 candidates
 963       FillTPCHistogramsNsigma(fListQAtpcV0,42,track,mult);
 964
 965     } // - End: track loop for pions from V0 -
 966
 967
 968     // - Begin: track loop for kaons from V0 -
 969     for(Int_t itrack = 0; itrack < fV0kaons->GetEntries(); itrack++){
 970       AliVTrack *track=(AliVTrack*)fV0kaons->At(itrack);
 971
 972       //
 973       //basic track cuts
 974       //
 975       ULong_t status=track->GetStatus();
 976       // not that nice. status bits not in virtual interface
 977       // TPC refit + ITS refit + TPC pid
 978       if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 979           !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
 980
 981       // The TPC pid cut removes the light nuclei (>5 sigma from proton line)
 982       //||        !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  )
 983       Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 984       Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 985       if (track->GetTPCNclsF()>0) {
 986         ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 987       }
 988
 989       if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 990
 991       eta=track->Eta();
 992       if ( TMath::Abs(eta)>0.9 ) continue;
 993
 994       // fill histograms for V0 candidates
 995       FillTPCHistogramsNsigma(fListQAtpcV0,43,track,mult);
 996
 997     } // - End: track loop for kaons from V0 -
 998
 999
1000     // - Begin: track loop for protons from V0 -
1001     for(Int_t itrack = 0; itrack < fV0protons->GetEntries(); itrack++){
1002       AliVTrack *track=(AliVTrack*)fV0protons->At(itrack);
1003
1004       //
1005       //basic track cuts
1006       //
1007       ULong_t status=track->GetStatus();
1008       // not that nice. status bits not in virtual interface
1009       // TPC refit + ITS refit + TPC pid
1010       if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
1011           !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
1012
1013       // The TPC pid cut removes the light nuclei (>5 sigma from proton line)
1014       //||        !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  )
1015       Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
1016       Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1017       if (track->GetTPCNclsF()>0) {
1018         ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
1019       }
1020
1021       if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
1022
1023       eta=track->Eta();
1024       if ( TMath::Abs(eta)>0.9 ) continue;
1025
1026       // fill histograms for V0 candidates
1027       FillTPCHistogramsNsigma(fListQAtpcV0,44,track,mult);
1028
1029     } // - End: track loop for protons from V0 -
1030
1031   } // -- End: track loops for V0 candidates --
1032
1033
1034   // Multiplicity distribution
1035   TH1 *hm=(TH1*)fListQAtpc->At(scCounter);
1036   if (hm) {
1037     hm->Fill(mult);
1038   }
1039
1040   //printf("\nAccepted multiplicity = %i \n --- END of event --- \n",nacc);
1041
1042 }
1043
1044 //______________________________________________________________________________
1045 void AliAnalysisTaskPIDqa::FillTRDqa()
1046 {
1047   //
1048   // Fill PID qa histograms for the TRD
1049   //
1050   AliVEvent *event=InputEvent();
1051   Int_t ntracks = event->GetNumberOfTracks();
1052   for(Int_t itrack = 0; itrack <  ntracks; itrack++){
1053     AliVTrack *track = (AliVTrack *)event->GetTrack(itrack);
1054
1055     //
1056     //basic track cuts
1057     //
1058     ULong_t status=track->GetStatus();
1059     // not that nice. status bits not in virtual interface
1060     // TPC refit + ITS refit + TPC pid + TRD out
1061     if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
1062         !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
1063 //         !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  ) || //removes light nuclei. So it is out for the moment
1064         !( (status & AliVTrack::kTRDout  ) == AliVTrack::kTRDout  )) continue;
1065
1066     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
1067     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1068     if (track->GetTPCNclsF()>0) {
1069       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
1070     }
1071
1072     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
1073
1074     Double_t likelihoods[AliPID::kSPECIES];
1075     if(fPIDResponse->ComputeTRDProbability(track, AliPID::kSPECIES, likelihoods) != AliPIDResponse::kDetPidOk) continue;
1076     Int_t ntracklets = 0;
1077     Double_t momentum = -1.;
1078     for(Int_t itl = 0; itl < 6; itl++) {
1079       if(track->GetTRDmomentum(itl) > 0.) {
1080         ntracklets++;
1081         if(momentum < 0) momentum = track->GetTRDmomentum(itl);
1082       }
1083     }
1084
1085     for(Int_t ispecie = 0; ispecie < AliPID::kSPECIES; ispecie++){
1086       TH2F *hLike = (TH2F *)fListQAtrd->At(ntracklets*AliPID::kSPECIES+ispecie);
1087       if (hLike) hLike->Fill(momentum,likelihoods[ispecie]);
1088     }
1089
1090     //=== nSigma and signal ===
1091     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1092       TH2 *h=(TH2*)fListQAtrdNsig->At(ispecie);
1093       TH2 *hTPCTOF=(TH2*)fListQAtrdNsigTPCTOF->At(ispecie);
1094       if (!h || !hTPCTOF) continue;
1095       Float_t nSigmaTPC=fPIDResponse->NumberOfSigmas(AliPIDResponse::kTPC, track, (AliPID::EParticleType)ispecie);
1096       Float_t nSigmaTRD=fPIDResponse->NumberOfSigmas(AliPIDResponse::kTRD, track, (AliPID::EParticleType)ispecie);
1097       Float_t nSigmaTOF=fPIDResponse->NumberOfSigmas(AliPIDResponse::kTOF, track, (AliPID::EParticleType)ispecie);
1098       h->Fill(momentum,nSigmaTRD);
1099
1100       if (TMath::Abs(nSigmaTPC)<3 && TMath::Abs(nSigmaTOF)<3) {
1101         hTPCTOF->Fill(momentum,nSigmaTRD);
1102       }
1103     }
1104
1105     TH2 *h=(TH2*)fListQAtrdNsig->Last();
1106
1107     if (h) {
1108       Double_t sig=track->GetTRDsignal();
1109       h->Fill(momentum,sig);
1110     }
1111
1112   }
1113 }
1114
1115 //______________________________________________________________________________
1116 void AliAnalysisTaskPIDqa::FillTOFqa()
1117 {
1118   //
1119   // Fill TOF information
1120   //
1121   AliVEvent *event=InputEvent();
1122
1123   Int_t ntracks=event->GetNumberOfTracks();
1124   Int_t tracksAtTof = 0;
1125   for(Int_t itrack = 0; itrack < ntracks; itrack++){
1126     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
1127
1128     //
1129     //basic track cuts
1130     //
1131     ULong_t status=track->GetStatus();
1132     // TPC refit + ITS refit +
1133     // TOF out + kTIME
1134     // kTIME
1135     // (we don't use kTOFmismatch because it depends on TPC and kTOFpid because it prevents light nuclei
1136     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
1137         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
1138         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||
1139         //        !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||
1140         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;
1141
1142     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
1143     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1144     if (track->GetTPCNclsF()>0) {
1145       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
1146     }
1147
1148     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
1149
1150     tracksAtTof++;
1151
1152     Double_t mom=track->P();
1153
1154     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1155       TH2 *h=(TH2*)fListQAtof->At(ispecie);
1156       if (!h) continue;
1157       Double_t nSigma=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
1158       h->Fill(mom,nSigma);
1159     }
1160
1161     TH2 *h=(TH2*)fListQAtof->FindObject("hSigP_TOF");
1162     if (h) {
1163       Double_t sig=track->GetTOFsignal()/1000.;
1164       h->Fill(mom,sig);
1165     }
1166
1167     Int_t mask = fPIDResponse->GetTOFResponse().GetStartTimeMask(mom);
1168     ((TH1F*)fListQAtof->FindObject("hStartTimeMask_TOF"))->Fill((Double_t)(mask+0.5));
1169
1170     if (mom >= 0.75 && mom <= 1.25 ) {
1171       Double_t nsigma= fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)AliPID::kPion);
1172       if (mask == 0) {
1173         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-Fill"))->Fill(nsigma);
1174       } else if (mask == 1) {
1175         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-TOF"))->Fill(nsigma);
1176       } else if ( (mask == 2) || (mask == 4) || (mask == 6) ) {
1177         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-T0"))->Fill(nsigma);
1178       } else {
1179         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-Best"))->Fill(nsigma);
1180       }
1181       if (mask & 0x1) { //at least TOF-T0 present
1182         Double_t delta=0;
1183         (void)fPIDResponse->GetSignalDelta((AliPIDResponse::EDetector)AliPIDResponse::kTOF,track,(AliPID::EParticleType)AliPID::kPion,delta);
1184         ((TH1F*)fListQAtof->FindObject("hDelta_TOF_Pion"))->Fill(delta);
1185       }
1186     }
1187
1188     Double_t res = (Double_t)fPIDResponse->GetTOFResponse().GetStartTimeRes(mom);
1189     ((TH1F*)fListQAtof->FindObject("hStartTimeRes_TOF"))->Fill(res);
1190
1191     Double_t startTimeT0 = event->GetT0TOF(0);
1192     if (startTimeT0 < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeAC_T0"))->Fill(startTimeT0);
1193     else {
1194       startTimeT0 = event->GetT0TOF(1);
1195       if (startTimeT0 < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeA_T0"))->Fill(startTimeT0);
1196       startTimeT0 = event->GetT0TOF(2);
1197       if (startTimeT0 < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeC_T0"))->Fill(startTimeT0);
1198     }
1199   }
1200   if (tracksAtTof > 0) {
1201     ((TH1F* )fListQAtof->FindObject("hnTracksAt_TOF"))->Fill(tracksAtTof);
1202     Int_t mask = fPIDResponse->GetTOFResponse().GetStartTimeMask(5.);
1203     if (mask & 0x1) ((TH1F*)fListQAtof->FindObject("hT0MakerEff"))->Fill(tracksAtTof);
1204   }
1205 }
1206
1207 //______________________________________________________________________________
1208 void AliAnalysisTaskPIDqa::FillT0qa()
1209 {
1210   //
1211   // Fill TOF information
1212   //
1213   AliVEvent *event=InputEvent();
1214
1215   Int_t ntracks=event->GetNumberOfTracks();
1216
1217   Int_t tracksAtT0 = 0;
1218
1219   for(Int_t itrack = 0; itrack < ntracks; itrack++){
1220     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
1221
1222     //
1223     //basic track cuts
1224     //
1225     ULong_t status=track->GetStatus();
1226     // TPC refit + ITS refit +
1227     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
1228         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
1229     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
1230     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1231     if (track->GetTPCNclsF()>0) {
1232       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
1233     }
1234     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
1235
1236     tracksAtT0++;
1237   }
1238
1239   Bool_t t0A = kFALSE;
1240   Bool_t t0C = kFALSE;
1241   Bool_t t0And = kFALSE;
1242   Double_t startTimeT0 = event->GetT0TOF(0);     // AND
1243   if (startTimeT0 < 90000) {
1244     t0And = kTRUE;
1245     ((TH1F*)fListQAt0->FindObject("hStartTimeAC_T0"))->Fill(startTimeT0);
1246     }
1247   startTimeT0 = event->GetT0TOF(1);             // T0A
1248   if (startTimeT0 < 90000) {
1249     t0A = kTRUE;
1250     ((TH1F*)fListQAt0->FindObject("hStartTimeA_T0"))->Fill(startTimeT0);
1251
1252   }
1253   startTimeT0 = event->GetT0TOF(2);             // T0C
1254   if (startTimeT0 < 90000) {
1255     t0C = kTRUE;
1256     ((TH1F*)fListQAt0->FindObject("hStartTimeC_T0"))->Fill(startTimeT0);
1257   }
1258
1259   ((TH1F* )fListQAt0->FindObject("hnTracksAt_T0"))->Fill(tracksAtT0);
1260   if (t0A) ((TH1F*)fListQAt0->FindObject("hT0AEff"))->Fill(tracksAtT0);
1261   if (t0C) ((TH1F*)fListQAt0->FindObject("hT0CEff"))->Fill(tracksAtT0);
1262   if (t0And) ((TH1F*)fListQAt0->FindObject("hT0AndEff"))->Fill(tracksAtT0);
1263   if (t0A || t0C) ((TH1F*)fListQAt0->FindObject("hT0OrEff"))->Fill(tracksAtT0);
1264 }
1265
1266
1267 //______________________________________________________________________________
1268 void AliAnalysisTaskPIDqa::FillEMCALqa()
1269 {
1270   //
1271   // Fill PID qa histograms for the EMCAL
1272   //
1273
1274   AliVEvent *event=InputEvent();
1275
1276   Int_t ntracks=event->GetNumberOfTracks();
1277   for(Int_t itrack = 0; itrack < ntracks; itrack++){
1278     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
1279
1280     //
1281     //basic track cuts
1282     //
1283     ULong_t status=track->GetStatus();
1284     // not that nice. status bits not in virtual interface
1285     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
1286
1287     Double_t pt=track->Pt();
1288
1289     //EMCAL nSigma (only for electrons at the moment)
1290     TH2 *h=(TH2*)fListQAemcal->At(0);
1291     if (!h) continue;
1292     Double_t nSigma=fPIDResponse->NumberOfSigmasEMCAL(track, (AliPID::EParticleType)0);
1293     h->Fill(pt,nSigma);
1294
1295   }
1296
1297    //EMCAL signal (E/p vs. pT) for electrons from V0
1298   for(Int_t itrack = 0; itrack < fV0electrons->GetEntries(); itrack++){
1299     AliVTrack *track=(AliVTrack*)fV0electrons->At(itrack);
1300
1301     //
1302     //basic track cuts
1303     //
1304     ULong_t status=track->GetStatus();
1305     // not that nice. status bits not in virtual interface
1306     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
1307
1308     Double_t pt=track->Pt();
1309
1310     TH2 *h=(TH2*)fListQAemcal->At(1);
1311     if (h) {
1312
1313       Int_t nMatchClus = track->GetEMCALcluster();
1314       Double_t mom     = track->P();
1315       Double_t eop     = -1.;
1316
1317       if(nMatchClus > -1){
1318
1319         AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
1320
1321         if(matchedClus){
1322
1323           // matched cluster is EMCAL
1324           if(matchedClus->IsEMCAL()){
1325
1326             Double_t fClsE       = matchedClus->E();
1327             eop                  = fClsE/mom;
1328
1329             h->Fill(pt,eop);
1330
1331           }
1332         }
1333       }
1334     }
1335   }
1336
1337    //EMCAL signal (E/p vs. pT) for pions from V0
1338   for(Int_t itrack = 0; itrack < fV0pions->GetEntries(); itrack++){
1339     AliVTrack *track=(AliVTrack*)fV0pions->At(itrack);
1340
1341     //
1342     //basic track cuts
1343     //
1344     ULong_t status=track->GetStatus();
1345     // not that nice. status bits not in virtual interface
1346     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
1347
1348     Double_t pt=track->Pt();
1349
1350     TH2 *h=(TH2*)fListQAemcal->At(2);
1351     if (h) {
1352
1353       Int_t nMatchClus = track->GetEMCALcluster();
1354       Double_t mom     = track->P();
1355       Double_t eop     = -1.;
1356
1357       if(nMatchClus > -1){
1358
1359         AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
1360
1361         if(matchedClus){
1362
1363           // matched cluster is EMCAL
1364           if(matchedClus->IsEMCAL()){
1365
1366             Double_t fClsE       = matchedClus->E();
1367             eop                  = fClsE/mom;
1368
1369             h->Fill(pt,eop);
1370
1371           }
1372         }
1373       }
1374     }
1375   }
1376
1377    //EMCAL signal (E/p vs. pT) for protons from V0
1378   for(Int_t itrack = 0; itrack < fV0protons->GetEntries(); itrack++){
1379     AliVTrack *track=(AliVTrack*)fV0protons->At(itrack);
1380
1381     //
1382     //basic track cuts
1383     //
1384     ULong_t status=track->GetStatus();
1385     // not that nice. status bits not in virtual interface
1386     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
1387
1388     Double_t pt=track->Pt();
1389
1390     TH2 *hP=(TH2*)fListQAemcal->At(3);
1391     TH2 *hAP=(TH2*)fListQAemcal->At(4);
1392     if (hP && hAP) {
1393
1394       Int_t nMatchClus = track->GetEMCALcluster();
1395       Double_t mom     = track->P();
1396       Int_t charge     = track->Charge();
1397       Double_t eop     = -1.;
1398
1399       if(nMatchClus > -1){
1400
1401         AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
1402
1403         if(matchedClus){
1404
1405           // matched cluster is EMCAL
1406           if(matchedClus->IsEMCAL()){
1407
1408             Double_t fClsE       = matchedClus->E();
1409             eop                  = fClsE/mom;
1410
1411             if(charge > 0)      hP->Fill(pt,eop);
1412             else if(charge < 0) hAP->Fill(pt,eop);
1413
1414           }
1415         }
1416       }
1417     }
1418   }
1419
1420 }
1421
1422
1423 //______________________________________________________________________________
1424 void AliAnalysisTaskPIDqa::FillHMPIDqa()
1425 {
1426   //
1427   // Fill PID qa histograms for the HMPID
1428   //
1429
1430   AliVEvent *event=InputEvent();
1431
1432   Int_t ntracks=event->GetNumberOfTracks();
1433   for(Int_t itrack = 0; itrack < ntracks; itrack++){
1434     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
1435
1436     //
1437     //basic track cuts
1438     //
1439     const ULong_t status=track->GetStatus();
1440     // not that nice. status bits not in virtual interface
1441     // TPC refit + ITS refit +
1442     // TOF out + TOFpid +
1443     // kTIME
1444     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
1445         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
1446
1447     const Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
1448     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1449     if (track->GetTPCNclsF()>0) {
1450       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
1451     }
1452
1453     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
1454
1455     const Double_t mom = track->P();
1456     const Double_t ckovAngle = track->GetHMPIDsignal();
1457
1458     Int_t nhists=0;
1459     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1460       if (ispecie==AliPID::kElectron || ispecie==AliPID::kMuon) continue;
1461       TH2 *h=(TH2*)fListQAhmpid->At(nhists);
1462       if (!h) {++nhists; continue;}
1463       const Double_t nSigma=fPIDResponse->NumberOfSigmasHMPID(track, (AliPID::EParticleType)ispecie);
1464       h->Fill(mom,nSigma);
1465       ++nhists;
1466     }
1467
1468     TH1F *hThetavsMom = (TH1F*)fListQAhmpid->At(AliPID::kSPECIESC);
1469
1470     if (hThetavsMom) hThetavsMom->Fill(mom,ckovAngle);
1471
1472   }
1473 }
1474 //______________________________________________________________________________
1475 void AliAnalysisTaskPIDqa::FillTOFHMPIDqa()
1476 {
1477   //
1478   // Fill PID qa histograms for the HMPID
1479   //
1480
1481   AliVEvent *event=InputEvent();
1482
1483   Int_t ntracks=event->GetNumberOfTracks();
1484   for(Int_t itrack = 0; itrack < ntracks; itrack++){
1485     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
1486
1487     //
1488     //basic track cuts
1489     //
1490     ULong_t status=track->GetStatus();
1491     // not that nice. status bits not in virtual interface
1492     // TPC refit + ITS refit +
1493     // TOF out + TOFpid +
1494     // kTIME
1495     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
1496         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
1497         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||
1498         !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||
1499         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;
1500
1501     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
1502     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1503     if (track->GetTPCNclsF()>0) {
1504       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
1505     }
1506
1507     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
1508
1509     Double_t mom = track->P();
1510     Double_t ckovAngle = track->GetHMPIDsignal();
1511
1512     Double_t nSigmaTOF[3];
1513     TH1F *h[3];
1514
1515     for (Int_t ispecie=2; ispecie<5; ++ispecie){
1516       //TOF nSigma
1517       nSigmaTOF[ispecie-2]=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
1518       h[ispecie-2] = (TH1F*)fListQAtofhmpid->At(ispecie-2);}
1519
1520     if(TMath::Abs(nSigmaTOF[0])<2)                                                              h[0]->Fill(mom,ckovAngle);
1521
1522     if(TMath::Abs(nSigmaTOF[1])<2 && TMath::Abs(nSigmaTOF[0])>3)                                h[1]->Fill(mom,ckovAngle);
1523
1524     if(TMath::Abs(nSigmaTOF[2])<2 && TMath::Abs(nSigmaTOF[1])>3 && TMath::Abs(nSigmaTOF[0])>3)  h[2]->Fill(mom,ckovAngle);
1525
1526   }
1527
1528 }
1529
1530 //______________________________________________________________________________
1531 void AliAnalysisTaskPIDqa::FillTPCTOFqa()
1532 {
1533   //
1534   // Fill PID qa histograms for the TOF
1535   //   Here also the TPC histograms after TOF selection are filled
1536   //
1537
1538   AliVEvent *event=InputEvent();
1539
1540   Int_t ntracks=event->GetNumberOfTracks();
1541   for(Int_t itrack = 0; itrack < ntracks; itrack++){
1542     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
1543
1544     //
1545     //basic track cuts
1546     //
1547     ULong_t status=track->GetStatus();
1548     // not that nice. status bits not in virtual interface
1549     // TPC refit + ITS refit +
1550     // TOF out + TOFpid +
1551     // kTIME
1552     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
1553         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
1554 //         !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid ) || //removes light nuclei, so it is out for the moment
1555         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||
1556         !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||
1557         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;
1558
1559     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
1560     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1561     if (track->GetTPCNclsF()>0) {
1562       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
1563     }
1564
1565     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
1566
1567
1568     Double_t mom=track->P();
1569     Double_t momTPC=track->GetTPCmomentum();
1570
1571     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1572       //TOF nSigma
1573       Double_t nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
1574       Double_t nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
1575
1576       //TPC after TOF cut
1577       TH2 *h=(TH2*)fListQAtpctof->At(ispecie);
1578       if (h && TMath::Abs(nSigmaTOF)<3.) h->Fill(momTPC,nSigmaTPC);
1579
1580       //TOF after TPC cut
1581       h=(TH2*)fListQAtpctof->At(ispecie+AliPID::kSPECIESC);
1582       if (h && TMath::Abs(nSigmaTPC)<3.) h->Fill(mom,nSigmaTOF);
1583
1584       //EMCAL after TOF and TPC cut
1585       h=(TH2*)fListQAtpctof->At(ispecie+2*AliPID::kSPECIESC);
1586       if (h && TMath::Abs(nSigmaTOF)<3. && TMath::Abs(nSigmaTPC)<3. ){
1587
1588         Int_t nMatchClus = track->GetEMCALcluster();
1589         Double_t pt      = track->Pt();
1590         Double_t eop     = -1.;
1591
1592         if(nMatchClus > -1){
1593
1594           AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
1595
1596           if(matchedClus){
1597
1598             // matched cluster is EMCAL
1599             if(matchedClus->IsEMCAL()){
1600
1601               Double_t fClsE       = matchedClus->E();
1602               eop                  = fClsE/mom;
1603
1604               h->Fill(pt,eop);
1605
1606
1607             }
1608           }
1609         }
1610       }
1611     }
1612   }
1613 }
1614
1615 //_____________________________________________________________________________
1616 void AliAnalysisTaskPIDqa::FillQAinfo()
1617 {
1618   //
1619   // Fill the QA information
1620   //
1621
1622
1623   //TPC QA info
1624   TObjArray *arrTPC=static_cast<TObjArray*>(fListQAinfo->At(0));
1625   if (fPIDResponse && arrTPC){
1626     AliTPCPIDResponse &tpcResp=fPIDResponse->GetTPCResponse();
1627     // fill spline names
1628     if (!arrTPC->UncheckedAt(0)){
1629
1630       TObjArray *arrTPCsplineNames=new TObjArray(AliPID::kSPECIESC);
1631       arrTPCsplineNames->SetOwner();
1632       arrTPCsplineNames->SetName("TPC_spline_names");
1633       arrTPC->AddAt(arrTPCsplineNames,0);
1634
1635       for (Int_t iresp=0; iresp<AliPID::kSPECIESC; ++iresp){
1636         const TObject *o=tpcResp.GetResponseFunction((AliPID::EParticleType)iresp);
1637         if (!o) continue;
1638         arrTPCsplineNames->Add(new TObjString(Form("%02d: %s",iresp, o->GetName())));
1639       }
1640     }
1641
1642     // tpc response config
1643     if (!arrTPC->UncheckedAt(1)){
1644
1645       TObjArray *arrTPCconfigInfo=new TObjArray;
1646       arrTPCconfigInfo->SetOwner();
1647       arrTPCconfigInfo->SetName("TPC_config_info");
1648       arrTPC->AddAt(arrTPCconfigInfo,1);
1649
1650       TObjString *ostr=0x0;
1651       ostr=new TObjString;
1652       ostr->String().Form("Eta Corr map: %s", tpcResp.GetEtaCorrMap()?tpcResp.GetEtaCorrMap()->GetName():"none");
1653       arrTPCconfigInfo->Add(ostr);
1654
1655       ostr=new TObjString;
1656       ostr->String().Form("Sigma Par map: %s", tpcResp.GetSigmaPar1Map()?tpcResp.GetSigmaPar1Map()->GetName():"none");
1657       arrTPCconfigInfo->Add(ostr);
1658
1659       ostr=new TObjString;
1660       ostr->String().Form("MIP: %.2f", tpcResp.GetMIP());
1661       arrTPCconfigInfo->Add(ostr);
1662
1663       ostr=new TObjString;
1664       ostr->String().Form("Res: Def %.3g (%.3g) : AllHigh %.3g (%.3g) : OROC high %.3g (%.3g)",
1665                           tpcResp.GetRes0(AliTPCPIDResponse::kDefault), tpcResp.GetResN2(AliTPCPIDResponse::kDefault),
1666                           tpcResp.GetRes0(AliTPCPIDResponse::kALLhigh), tpcResp.GetResN2(AliTPCPIDResponse::kALLhigh),
1667                           tpcResp.GetRes0(AliTPCPIDResponse::kOROChigh), tpcResp.GetResN2(AliTPCPIDResponse::kOROChigh)
1668                          );
1669       arrTPCconfigInfo->Add(ostr);
1670     }
1671   }
1672 }
1673
1674 //______________________________________________________________________________
1675 void AliAnalysisTaskPIDqa::SetupITSqa()
1676 {
1677   //
1678   // Create the ITS qa objects
1679   //
1680
1681   TVectorD *vX=MakeLogBinning(200,.1,30);
1682
1683   //ITS+TPC tracks
1684   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1685     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_ITS_%s",AliPID::ParticleName(ispecie)),
1686                               Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1687                               vX->GetNrows()-1,vX->GetMatrixArray(),
1688                               200,-10,10);
1689     fListQAits->Add(hNsigmaP);
1690   }
1691   TH2F *hSig = new TH2F("hSigP_ITS",
1692                         "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
1693                         vX->GetNrows()-1,vX->GetMatrixArray(),
1694                         300,0,300);
1695   fListQAits->Add(hSig);
1696
1697   //ITS Standalone tracks
1698   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1699     TH2F *hNsigmaPSA = new TH2F(Form("hNsigmaP_ITSSA_%s",AliPID::ParticleName(ispecie)),
1700                                 Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1701                                 vX->GetNrows()-1,vX->GetMatrixArray(),
1702                                 200,-10,10);
1703     fListQAitsSA->Add(hNsigmaPSA);
1704   }
1705   TH2F *hSigSA = new TH2F("hSigP_ITSSA",
1706                           "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
1707                           vX->GetNrows()-1,vX->GetMatrixArray(),
1708                           300,0,300);
1709   fListQAitsSA->Add(hSigSA);
1710
1711   //ITS Pure Standalone tracks
1712   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1713     TH2F *hNsigmaPPureSA = new TH2F(Form("hNsigmaP_ITSPureSA_%s",AliPID::ParticleName(ispecie)),
1714                                     Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1715                                     vX->GetNrows()-1,vX->GetMatrixArray(),
1716                                     200,-10,10);
1717     fListQAitsPureSA->Add(hNsigmaPPureSA);
1718   }
1719   TH2F *hSigPureSA = new TH2F("hSigP_ITSPureSA",
1720                               "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
1721                               vX->GetNrows()-1,vX->GetMatrixArray(),
1722                               300,0,300);
1723   fListQAitsPureSA->Add(hSigPureSA);
1724
1725   delete vX;
1726 }
1727
1728 //_____________________________________________________________________________
1729 void AliAnalysisTaskPIDqa::AddTPCHistogramsSignal(TList *sublist, const char *scenario)
1730 {
1731   //
1732   // Create the TPC qa objects: create histograms for the TPC signal for different settings
1733   //
1734
1735   TVectorD *vX=MakeLogBinning(200,.1,30);
1736   Int_t nBinsMult = 38;
1737   Double_t xBinsMult[39] = {0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100,
1738                            120, 140, 160, 180, 200,
1739                            300, 400, 500, 600, 700, 800, 900, 1000,
1740                            1200, 1400, 1600, 1800, 2000,
1741                            2200, 2400, 2600, 2800, 3000,
1742                            3200, 3400, 3600, 3800, 4000
1743                            };
1744   const Int_t binsEta=110;
1745   Float_t etaMin=-1.1;
1746   Float_t etaMax=1.1;
1747
1748   char signal[4][12]={"std","IROC","OROCmedium","OROClong"};
1749
1750
1751   // TPC signal vs. p for all particles (standard, IROC, OROCmedium, OROClong)
1752   for (Int_t iSig=0; iSig<4; iSig++) {
1753     TH2F *hSigP = new TH2F(Form("hSigP_TPC_%s_%s",signal[iSig],scenario),
1754                              Form("TPC_%s signal (%s) vs. p;p [GeV]; TPC signal [arb. units]",scenario,signal[iSig]),
1755                              vX->GetNrows()-1,vX->GetMatrixArray(),
1756                              300,0,300);
1757     sublist->Add(hSigP);
1758   }
1759
1760   // MIP pions: TPC signal vs. eta
1761   for (Int_t iSig=0; iSig<4; iSig++) {
1762     TH2F *hSigEtaMIPpi = new TH2F(Form("hSigEta_TPC_%s_%s_MIPpi",signal[iSig],scenario),
1763                                   Form("TPC_%s signal (%s) MIPpi vs. eta;#eta;TPC signal [arb. units]",scenario,signal[iSig]),
1764                                   binsEta,etaMin,etaMax,
1765                                   300,0,300);
1766     sublist->Add(hSigEtaMIPpi);
1767   }
1768
1769   // MIP pions: TPC signal vs. multiplicity
1770   for (Int_t iSig=0; iSig<4; iSig++) {
1771     TH2F *hSigMultMPIpi = new TH2F(Form("hSigMult_TPC_%s_%s_MIPpi",signal[iSig],scenario),
1772                                    Form("TPC_%s signal (%s) MIPpi vs. mult;multiplicity;TPC signal [arb. units]",scenario,signal[iSig]),
1773                                    nBinsMult,xBinsMult,
1774                                    300,0,300);
1775     sublist->Add(hSigMultMPIpi);
1776   }
1777
1778   // Electrons: TPC signal vs. eta
1779   for (Int_t iSig=0; iSig<4; iSig++) {
1780     TH2F *hSigEtaEle = new TH2F(Form("hSigEta_TPC_%s_%s_Ele",signal[iSig],scenario),
1781                                 Form("TPC_%s signal (%s) electrons vs. eta;#eta;TPC signal [arb. units]",scenario,signal[iSig]),
1782                                 binsEta,etaMin,etaMax,
1783                                 300,0,300);
1784     sublist->Add(hSigEtaEle);
1785   }
1786
1787   // Electrons: TPC signal vs. multiplicity
1788   for (Int_t iSig=0; iSig<4; iSig++) {
1789     TH2F *hSigMultEle = new TH2F(Form("hSigMult_TPC_%s_%s_Ele",signal[iSig],scenario),
1790                                  Form("TPC_%s signal (%s) electrons vs. mult;multiplicity;TPC signal [arb. units]",scenario,signal[iSig]),
1791                                  nBinsMult,xBinsMult,
1792                                  300,0,300);
1793     sublist->Add(hSigMultEle);
1794   }
1795
1796   delete vX;
1797
1798 }
1799
1800 //_____________________________________________________________________________
1801 void AliAnalysisTaskPIDqa::AddTPCHistogramsNsigma(TList *sublist, const char *scenario, Int_t scnumber)
1802 {
1803   //
1804   // Create the TPC qa objects: create histograms for TPC Nsigma for different settings
1805   //
1806
1807   TVectorD *vX=MakeLogBinning(200,.1,30.);
1808   Int_t nBinsMult = 38;
1809   Double_t xBinsMult[39] = {0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100,
1810                            120, 140, 160, 180, 200,
1811                            300, 400, 500, 600, 700, 800, 900, 1000,
1812                            1200, 1400, 1600, 1800, 2000,
1813                            2200, 2400, 2600, 2800, 3000,
1814                            3200, 3400, 3600, 3800, 4000
1815                            };
1816   const Int_t binsEta=110;
1817   Float_t etaMin=-1.1;
1818   Float_t etaMax=1.1;
1819
1820   Int_t nSpecies=0;
1821
1822   if (scnumber == 4) nSpecies=(Int_t)AliPID::kSPECIES;
1823   else nSpecies=(Int_t)AliPID::kSPECIESC;
1824
1825   // Nsigma vs. p for different particle species
1826   for (Int_t ispecie=0; ispecie<nSpecies; ++ispecie){
1827     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_%s_%s",scenario,AliPID::ParticleName(ispecie)),
1828                               Form("TPC_%s n#sigma %s vs. p;p [GeV]; n#sigma",scenario,AliPID::ParticleName(ispecie)),
1829                               vX->GetNrows()-1,vX->GetMatrixArray(),
1830                               200,-10,10);
1831     sublist->Add(hNsigmaP);
1832   }
1833
1834   // Nsigma vs. eta for different particle species (only for some scenarios)
1835   if ( scnumber == 1 || scnumber == 4 ) {
1836     for (Int_t ispecie=0; ispecie<nSpecies; ++ispecie){
1837       TH2F *hNsigmaEta = new TH2F(Form("hNsigmaEta_TPC_%s_%s",scenario,AliPID::ParticleName(ispecie)),
1838                               Form("TPC_%s n#sigma %s vs. eta;#eta; n#sigma",scenario,AliPID::ParticleName(ispecie)),
1839                               binsEta,etaMin,etaMax,
1840                               200,-10,10);
1841       sublist->Add(hNsigmaEta);
1842     }
1843   }
1844
1845   // Nsigma vs. multiplicity for different particle species (only for some scenarios)
1846   if ( scnumber == 1 || scnumber == 4 ) {
1847     for (Int_t ispecie=0; ispecie<nSpecies; ++ispecie){
1848       TH2F *hNsigmaMult = new TH2F(Form("hNsigmaMult_TPC_%s_%s",scenario,AliPID::ParticleName(ispecie)),
1849                               Form("TPC_%s n#sigma %s vs. mult;multiplicity; n#sigma",scenario,AliPID::ParticleName(ispecie)),
1850                               nBinsMult,xBinsMult,
1851                               200,-10,10);
1852       sublist->Add(hNsigmaMult);
1853     }
1854   }
1855
1856   // - Beginn: Adding histograms for MIP pions and electrons (only for some scenarios) -
1857   if ( scnumber == 0 || scnumber == 2 || scnumber == 3 ) {
1858
1859     // MIP pions: Nsigma vs. eta
1860     TH2F *hNsigmaEtaMIPpi = new TH2F(Form("hNsigmaEta_TPC_%s_MIPpi",scenario),
1861                               Form("TPC_%s n#sigma MIPpi vs. eta;#eta; n#sigma",scenario),
1862                               binsEta,etaMin,etaMax,
1863                               200,-10,10);
1864     sublist->Add(hNsigmaEtaMIPpi);
1865
1866     // MIP pions: Nsigma vs. multiplicity
1867     TH2F *hNsigmaMultMIPpi = new TH2F(Form("hNsigmaMult_TPC_%s_MIPpi",scenario),
1868                                Form("TPC_%s n#sigma MIPpi vs. mult;multiplicity; n#sigma",scenario),
1869                                nBinsMult,xBinsMult,
1870                                200,-10,10);
1871     sublist->Add(hNsigmaMultMIPpi);
1872
1873     // Electrons: Nsigma vs. eta
1874     TH2F *hNsigmaEtaEle = new TH2F(Form("hNsigmaEta_TPC_%s_Ele",scenario),
1875                               Form("TPC_%s n#sigma electrons vs. eta;#eta; n#sigma",scenario),
1876                               binsEta,etaMin,etaMax,
1877                               200,-10,10);
1878     sublist->Add(hNsigmaEtaEle);
1879
1880     // Electrons: Nsigma vs. multiplicity
1881     TH2F *hNsigmaMultEle = new TH2F(Form("hNsigmaMult_TPC_%s_Ele",scenario),
1882                                Form("TPC_%s n#sigma electrons vs. mult;multiplicity; n#sigma",scenario),
1883                                nBinsMult,xBinsMult,
1884                                200,-10,10);
1885     sublist->Add(hNsigmaMultEle);
1886   } // - End: Adding histograms for MIP pions and electrons
1887
1888   delete vX;
1889
1890 }
1891
1892 //______________________________________________________________________________
1893 void AliAnalysisTaskPIDqa::SetupTPCqa(Bool_t fillMC, Bool_t fill11h, Bool_t fillV0)
1894 {
1895   //
1896   // Create the TPC qa objects
1897   //
1898
1899   // Set up the multiplicity binning
1900   Int_t nBinsMult = 38;
1901   Double_t xBinsMult[39] = {0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100,
1902                            120, 140, 160, 180, 200,
1903                            300, 400, 500, 600, 700, 800, 900, 1000,
1904                            1200, 1400, 1600, 1800, 2000,
1905                            2200, 2400, 2600, 2800, 3000,
1906                            3200, 3400, 3600, 3800, 4000
1907                            };
1908
1909
1910   // Create TPC sublists for different scenarios
1911   // corresponding to available information,
1912   // e.g. MC or not, special settings for LHC11h
1913
1914   // basic/default scenario, used always
1915   fListQAtpcBasic=new TList;
1916   fListQAtpcBasic->SetOwner();
1917   fListQAtpcBasic->SetName("TPCBasic");
1918   fListQAtpc->Add(fListQAtpcBasic);
1919
1920   // MC truth scenario: use only MC truth identified particles
1921   // only available for MC
1922   if (fillMC == kTRUE) {
1923     fListQAtpcMCtruth=new TList;
1924     fListQAtpcMCtruth->SetOwner();
1925     fListQAtpcMCtruth->SetName("TPCMCtruth");
1926     fListQAtpc->Add(fListQAtpcMCtruth);
1927   }
1928
1929   // Hybrid and OROChigh scenarios,
1930   // special settings only available for PbPb LHC11h data
1931   if (fill11h == kTRUE) {
1932     fListQAtpcHybrid=new TList;
1933     fListQAtpcHybrid->SetOwner();
1934     fListQAtpcHybrid->SetName("TPCHybrid");
1935     fListQAtpc->Add(fListQAtpcHybrid);
1936
1937     fListQAtpcOROChigh=new TList;
1938     fListQAtpcOROChigh->SetOwner();
1939     fListQAtpcOROChigh->SetName("TPCOROChigh");
1940     fListQAtpc->Add(fListQAtpcOROChigh);
1941   }
1942
1943   // scenario only for V0s,
1944   // only available for ESDs
1945   if (fillV0 == kTRUE) {
1946     fListQAtpcV0=new TList;
1947     fListQAtpcV0->SetOwner();
1948     fListQAtpcV0->SetName("TPCV0");
1949     fListQAtpc->Add(fListQAtpcV0);
1950   }
1951
1952
1953   // the default ("basic") scenario
1954   AddTPCHistogramsNsigma(fListQAtpcBasic,"Basic",0);
1955   AddTPCHistogramsSignal(fListQAtpcBasic,"Basic");
1956
1957   // only MC truth identified particles
1958   if (fillMC) {
1959     AddTPCHistogramsNsigma(fListQAtpcMCtruth,"MCtruth",1);
1960   }
1961
1962   // the "hybrid" scenario (only for period LHC11h)
1963   if (fill11h) {
1964     AddTPCHistogramsNsigma(fListQAtpcHybrid,"Hybrid",2);
1965   }
1966
1967   // the "OROC high" scenario (only for period LHC11h)
1968   if (fill11h) {
1969     AddTPCHistogramsNsigma(fListQAtpcOROChigh,"OROChigh",3);
1970   }
1971
1972   // only for V0s
1973   if (fillV0) {
1974     AddTPCHistogramsNsigma(fListQAtpcV0,"V0",4);
1975   }
1976
1977
1978   // Multiplicity distribution --- as check
1979   TH1F *hMult = new TH1F("hMult_TPC",
1980                          "Multiplicity distribution;multiplicity;counts",
1981                          nBinsMult,xBinsMult);
1982   fListQAtpc->Add(hMult);
1983
1984 }
1985
1986 //______________________________________________________________________________
1987 void AliAnalysisTaskPIDqa::SetupTRDqa()
1988 {
1989   //
1990   // Create the TRD qa objects
1991   //
1992   TVectorD *vX=MakeLogBinning(200,.1,30);
1993   for(Int_t itl = 0; itl < 6; ++itl){
1994     for(Int_t ispecie = 0; ispecie < AliPID::kSPECIES; ispecie++){
1995       TH2F *hLikeP = new TH2F(Form("hLikeP_TRD_%dtls_%s", itl, AliPID::ParticleName(ispecie)),
1996                               Form("TRD Likelihood to be %s %s for tracks having %d %s; p (GeV/c); TRD %s Likelihood", ispecie == 0 ? "an" : "a", AliPID::ParticleName(ispecie), itl+1, itl == 0 ? "tracklet" : "tracklets", AliPID::ParticleName(ispecie)),
1997                               vX->GetNrows()-1, vX->GetMatrixArray(),
1998                               100, 0., 1.);
1999       fListQAtrd->Add(hLikeP);
2000     }
2001   }
2002
2003   // === nSigma Values and signal ===
2004   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
2005     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TRD_%s",AliPID::ParticleName(ispecie)),
2006                               Form("TRD n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
2007                               vX->GetNrows()-1,vX->GetMatrixArray(),
2008                               100,-10,10);
2009     fListQAtrdNsig->Add(hNsigmaP);
2010   }
2011
2012   TH2F *hSig = new TH2F("hSigP_TRD",
2013                         "TRD signal vs. p;p [GeV]; TRD signal [arb. units]",
2014                         vX->GetNrows()-1,vX->GetMatrixArray(),
2015                         100,0,100);
2016   fListQAtrdNsig->Add(hSig);
2017
2018   fListQAtrd->Add(fListQAtrdNsig);
2019
2020   // === Same after 3 sigma in TPC and TOF
2021   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
2022     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TRD_TPCTOF_%s",AliPID::ParticleName(ispecie)),
2023                               Form("TRD n#sigma %s vs. p after 3#sigma cut in TPC&TOF;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
2024                               vX->GetNrows()-1,vX->GetMatrixArray(),
2025                               100,-10,10);
2026     fListQAtrdNsigTPCTOF->Add(hNsigmaP);
2027   }
2028
2029   fListQAtrd->Add(fListQAtrdNsigTPCTOF);
2030
2031   delete vX;
2032 }
2033
2034 //______________________________________________________________________________
2035 void AliAnalysisTaskPIDqa::SetupTOFqa()
2036 {
2037   //
2038   // Create the TOF qa objects
2039   //
2040
2041   TVectorD *vX=MakeLogBinning(200,.1,30);
2042
2043   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
2044     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_%s",AliPID::ParticleName(ispecie)),
2045                               Form("TOF n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
2046                               vX->GetNrows()-1,vX->GetMatrixArray(),
2047                               200,-10,10);
2048     fListQAtof->Add(hNsigmaP);
2049   }
2050
2051   TH1F *hnSigT0Fill = new TH1F("hNsigma_TOF_Pion_T0-Fill","TOF n#sigma (Pion) T0-FILL [0.75-1.25. GeV/c]",200,-10,10);
2052   fListQAtof->Add(hnSigT0Fill);
2053   TH1F *hnSigT0T0 = new TH1F("hNsigma_TOF_Pion_T0-T0","TOF n#sigma (Pion) T0-T0 [0.75-1.25 GeV/c]",200,-10,10);
2054   fListQAtof->Add(hnSigT0T0);
2055   TH1F *hnSigT0TOF = new TH1F("hNsigma_TOF_Pion_T0-TOF","TOF n#sigma (Pion) T0-TOF [0.75-1.25 GeV/c]",200,-10,10);
2056   fListQAtof->Add(hnSigT0TOF);
2057   TH1F *hnSigT0Best = new TH1F("hNsigma_TOF_Pion_T0-Best","TOF n#sigma (Pion) T0-Best [0.75-1.25 GeV/c]",200,-10,10);
2058   fListQAtof->Add(hnSigT0Best);
2059   TH1F *hnDeltaPi = new TH1F("hDelta_TOF_Pion","DeltaT (Pion) [0.75-1.25 GeV/c]",50,-500,500);
2060   fListQAtof->Add(hnDeltaPi);
2061
2062   TH2F *hSig = new TH2F("hSigP_TOF",
2063                         "TOF signal vs. p;p [GeV]; TOF signal [ns]",
2064                         vX->GetNrows()-1,vX->GetMatrixArray(),
2065                         300,0,30);
2066
2067   delete vX;
2068
2069   fListQAtof->Add(hSig);
2070
2071   TH1F *hStartTimeMaskTOF = new TH1F("hStartTimeMask_TOF","StartTime mask",8,0,8);
2072   fListQAtof->Add(hStartTimeMaskTOF);
2073   TH1F *hStartTimeResTOF = new TH1F("hStartTimeRes_TOF","StartTime resolution [ps]",100,0,500);
2074   fListQAtof->Add(hStartTimeResTOF);
2075
2076   TH1F *hnTracksAtTOF = new TH1F("hnTracksAt_TOF","Matched tracks at TOF",100,0,100);
2077   fListQAtof->Add(hnTracksAtTOF);
2078   TH1F *hT0MakerEff = new TH1F("hT0MakerEff","Events with T0-TOF vs nTracks",100,0,100);
2079   fListQAtof->Add(hT0MakerEff);
2080
2081   // this in principle should stay on a T0 PID QA, but are just the data prepared for TOF use
2082   TH1F *hStartTimeAT0 = new TH1F("hStartTimeA_T0","StartTime from T0A [ps]",1000,-1000,1000);
2083   fListQAtof->Add(hStartTimeAT0);
2084   TH1F *hStartTimeCT0 = new TH1F("hStartTimeC_T0","StartTime from T0C [ps]",1000,-1000,1000);
2085   fListQAtof->Add(hStartTimeCT0);
2086   TH1F *hStartTimeACT0 = new TH1F("hStartTimeAC_T0","StartTime from T0AC [ps]",1000,-1000,1000);;
2087   fListQAtof->Add(hStartTimeACT0);
2088 }
2089
2090
2091 //______________________________________________________________________________
2092 void AliAnalysisTaskPIDqa::SetupT0qa()
2093 {
2094   //
2095   // Create the T0 qa objects
2096   //
2097
2098   // these are similar to plots inside TOFqa, but these are for all events
2099   TH1F *hStartTimeAT0 = new TH1F("hStartTimeA_T0","StartTime from T0A [ps]",1000,-1000,1000);
2100   fListQAt0->Add(hStartTimeAT0);
2101   TH1F *hStartTimeCT0 = new TH1F("hStartTimeC_T0","StartTime from T0C [ps]",1000,-1000,1000);
2102   fListQAt0->Add(hStartTimeCT0);
2103   TH1F *hStartTimeACT0 = new TH1F("hStartTimeAC_T0","StartTime from T0AC [ps]",1000,-1000,1000);;
2104   fListQAt0->Add(hStartTimeACT0);
2105
2106   TH1F *hnTracksAtT0 = new TH1F("hnTracksAt_T0","Tracks for events selected for T0",100,0,100);
2107   fListQAt0->Add(hnTracksAtT0);
2108   TH1F *hT0AEff = new TH1F("hT0AEff","Events with T0A vs nTracks",100,0,100);
2109   fListQAt0->Add(hT0AEff);
2110   TH1F *hT0CEff = new TH1F("hT0CEff","Events with T0C vs nTracks",100,0,100);
2111   fListQAt0->Add(hT0CEff);
2112   TH1F *hT0AndEff = new TH1F("hT0AndEff","Events with T0AC (AND) vs nTracks",100,0,100);
2113   fListQAt0->Add(hT0AndEff);
2114   TH1F *hT0OrEff = new TH1F("hT0OrEff","Events with T0AC (OR) vs nTracks",100,0,100);
2115   fListQAt0->Add(hT0OrEff);
2116
2117
2118 }
2119
2120 //______________________________________________________________________________
2121 void AliAnalysisTaskPIDqa::SetupEMCALqa()
2122 {
2123   //
2124   // Create the EMCAL qa objects
2125   //
2126
2127   TVectorD *vX=MakeLogBinning(200,.1,30);
2128
2129   TH2F *hNsigmaPt = new TH2F(Form("hNsigmaPt_EMCAL_%s",AliPID::ParticleName(0)),
2130                              Form("EMCAL n#sigma %s vs. p_{T};p_{T} [GeV]; n#sigma",AliPID::ParticleName(0)),
2131                              vX->GetNrows()-1,vX->GetMatrixArray(),
2132                              200,-10,10);
2133   fListQAemcal->Add(hNsigmaPt);
2134
2135
2136   TH2F *hSigPtEle = new TH2F("hSigPt_EMCAL_Ele",
2137                         "EMCAL signal (E/p) vs. p_{T} for electrons;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
2138                         vX->GetNrows()-1,vX->GetMatrixArray(),
2139                         200,0,2);
2140   fListQAemcal->Add(hSigPtEle);
2141
2142   TH2F *hSigPtPions = new TH2F("hSigPt_EMCAL_Pions",
2143                         "EMCAL signal (E/p) vs. p_{T} for pions;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
2144                         vX->GetNrows()-1,vX->GetMatrixArray(),
2145                         200,0,2);
2146   fListQAemcal->Add(hSigPtPions);
2147
2148   TH2F *hSigPtProtons = new TH2F("hSigPt_EMCAL_Protons",
2149                         "EMCAL signal (E/p) vs. p_{T} for protons;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
2150                         vX->GetNrows()-1,vX->GetMatrixArray(),
2151                         200,0,2);
2152   fListQAemcal->Add(hSigPtProtons);
2153
2154   TH2F *hSigPtAntiProtons = new TH2F("hSigPt_EMCAL_Antiprotons",
2155                         "EMCAL signal (E/p) vs. p_{T} for antiprotons;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
2156                         vX->GetNrows()-1,vX->GetMatrixArray(),
2157                         200,0,2);
2158   fListQAemcal->Add(hSigPtAntiProtons);
2159
2160   delete vX;
2161 }
2162
2163 //______________________________________________________________________________
2164 void AliAnalysisTaskPIDqa::SetupHMPIDqa()
2165 {
2166   //
2167   // Create the HMPID qa objects
2168   //
2169
2170   TVectorD *vX=MakeLogBinning(200,.1,30);
2171
2172   // nSigmas
2173   Int_t nhists=0;
2174   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
2175     if (ispecie==AliPID::kElectron || ispecie==AliPID::kMuon) continue;
2176     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_HMPID_%s",AliPID::ParticleName(ispecie)),
2177                               Form("HMPID n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
2178                               vX->GetNrows()-1,vX->GetMatrixArray(),
2179                               200,-10,10);
2180     fListQAhmpid->AddAt(hNsigmaP, nhists);
2181     ++nhists;
2182   }
2183
2184   // cherenkov angle
2185   TH2F *hCkovAnglevsMom   = new TH2F("hCkovAnglevsMom",  "Cherenkov angle vs momentum",
2186                                      vX->GetNrows()-1,vX->GetMatrixArray(),
2187                                      500,0,1);
2188   fListQAhmpid->AddAt(hCkovAnglevsMom,nhists);
2189
2190   delete vX;
2191 }
2192
2193 //______________________________________________________________________________
2194 void AliAnalysisTaskPIDqa::SetupTOFHMPIDqa()
2195 {
2196   //
2197   // Create the HMPID qa objects
2198   //
2199
2200   TH2F *hCkovAnglevsMomPion   = new TH2F("hCkovAnglevsMom_pion",  "Cherenkov angle vs momentum for pions",500,0,5.,500,0,1);
2201   fListQAtofhmpid->Add(hCkovAnglevsMomPion);
2202
2203   TH2F *hCkovAnglevsMomKaon   = new TH2F("hCkovAnglevsMom_kaon",  "Cherenkov angle vs momentum for kaons",500,0,5.,500,0,1);
2204   fListQAtofhmpid->Add(hCkovAnglevsMomKaon);
2205
2206   TH2F *hCkovAnglevsMomProton = new TH2F("hCkovAnglevsMom_proton","Cherenkov angle vs momentum for protons",500,0,5.,500,0,1);
2207   fListQAtofhmpid->Add(hCkovAnglevsMomProton);
2208
2209
2210 }
2211
2212 //______________________________________________________________________________
2213 void AliAnalysisTaskPIDqa::SetupTPCTOFqa()
2214 {
2215   //
2216   // Create the qa objects for TPC + TOF combination
2217   //
2218
2219   TVectorD *vX=MakeLogBinning(200,.1,30);
2220
2221   //TPC signals after TOF cut
2222   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
2223     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_TOF_%s",AliPID::ParticleName(ispecie)),
2224                               Form("TPC n#sigma %s vs. p (after TOF 3#sigma cut);p_{TPC} [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
2225                               vX->GetNrows()-1,vX->GetMatrixArray(),
2226                               200,-10,10);
2227     fListQAtpctof->Add(hNsigmaP);
2228   }
2229
2230   //TOF signals after TPC cut
2231   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
2232     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_TPC_%s",AliPID::ParticleName(ispecie)),
2233                               Form("TOF n#sigma %s vs. p (after TPC n#sigma cut);p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
2234                               vX->GetNrows()-1,vX->GetMatrixArray(),
2235                               200,-10,10);
2236     fListQAtpctof->Add(hNsigmaP);
2237   }
2238
2239   //EMCAL signal after TOF and TPC cut
2240   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
2241     TH2F *heopPt = new TH2F(Form("heopPt_TOF_TPC_%s",AliPID::ParticleName(ispecie)),
2242                             Form("EMCAL signal (E/p) %s vs. p_{T};p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",AliPID::ParticleName(ispecie)),
2243                             vX->GetNrows()-1,vX->GetMatrixArray(),
2244                             200,0,2);
2245     fListQAtpctof->Add(heopPt);
2246   }
2247
2248   delete vX;
2249 }
2250 //______________________________________________________________________________
2251 void AliAnalysisTaskPIDqa::SetupV0qa()
2252 {
2253   //
2254   // Create the qa objects for V0 Kine cuts
2255   //
2256
2257   TH2F *hArmenteros  = new TH2F("hArmenteros",  "Armenteros plot",200,-1.,1.,200,0.,0.4);
2258   fListQAV0->Add(hArmenteros);
2259
2260 }
2261
2262 //_____________________________________________________________________________
2263 void AliAnalysisTaskPIDqa::SetupQAinfo(){
2264   //
2265   // Setup the info of QA objects
2266   //
2267
2268   TObjArray *arr=new TObjArray;
2269   arr->SetName("TPC_info");
2270   fListQAinfo->Add(arr);
2271 }
2272
2273 //______________________________________________________________________________
2274 TVectorD* AliAnalysisTaskPIDqa::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
2275 {
2276   //
2277   // Make logarithmic binning
2278   // the user has to delete the array afterwards!!!
2279   //
2280
2281   //check limits
2282   if (xmin<1e-20 || xmax<1e-20){
2283     AliError("For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
2284     return MakeLinBinning(nbinsX, xmin, xmax);
2285   }
2286   if (xmax<xmin){
2287     Double_t tmp=xmin;
2288     xmin=xmax;
2289     xmax=tmp;
2290   }
2291   TVectorD *binLim=new TVectorD(nbinsX+1);
2292   Double_t first=xmin;
2293   Double_t last=xmax;
2294   Double_t expMax=TMath::Log(last/first);
2295   for (Int_t i=0; i<nbinsX+1; ++i){
2296     (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
2297   }
2298   return binLim;
2299 }
2300
2301 //______________________________________________________________________________
2302 TVectorD* AliAnalysisTaskPIDqa::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
2303 {
2304   //
2305   // Make linear binning
2306   // the user has to delete the array afterwards!!!
2307   //
2308   if (xmax<xmin){
2309     Double_t tmp=xmin;
2310     xmin=xmax;
2311     xmax=tmp;
2312   }
2313   TVectorD *binLim=new TVectorD(nbinsX+1);
2314   Double_t first=xmin;
2315   Double_t last=xmax;
2316   Double_t binWidth=(last-first)/nbinsX;
2317   for (Int_t i=0; i<nbinsX+1; ++i){
2318     (*binLim)[i]=first+binWidth*(Double_t)i;
2319   }
2320   return binLim;
2321 }
2322
2323 //_____________________________________________________________________________
2324 TVectorD* AliAnalysisTaskPIDqa::MakeArbitraryBinning(const char* bins)
2325 {
2326   //
2327   // Make arbitrary binning, bins separated by a ','
2328   //
2329   TString limits(bins);
2330   if (limits.IsNull()){
2331     AliError("Bin Limit string is empty, cannot add the variable");
2332     return 0x0;
2333   }
2334
2335   TObjArray *arr=limits.Tokenize(",");
2336   Int_t nLimits=arr->GetEntries();
2337   if (nLimits<2){
2338     AliError("Need at leas 2 bin limits, cannot add the variable");
2339     delete arr;
2340     return 0x0;
2341   }
2342
2343   TVectorD *binLimits=new TVectorD(nLimits);
2344   for (Int_t iLim=0; iLim<nLimits; ++iLim){
2345     (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
2346   }
2347
2348   delete arr;
2349   return binLimits;
2350 }
2351