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     AliAODHeader * header=dynamic_cast<AliAODHeader*>(fAODevent->GetHeader());
 838     if(!header) AliFatal("Not a standard AOD");
 839     mult=header->GetTPConlyRefMultiplicity();
 840   }
 841
 842   /*if (mult < 0) {
 843     printf("Reference multiplicity not available \n");
 844     //return;
 845   }*/
 846
 847   //printf("The multiplicity is = %i ",mult);
 848
 849
 850   // -- Begin: track loop --
 851   Int_t ntracks=event->GetNumberOfTracks();
 852   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 853     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 854
 855     //
 856     //basic track cuts
 857     //
 858     ULong_t status=track->GetStatus();
 859     // not that nice. status bits not in virtual interface
 860     // TPC refit + ITS refit + TPC pid
 861     if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 862         !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
 863
 864     // The TPC pid cut removes the light nuclei (>5 sigma from proton line)
 865     //||        !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  )
 866     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 867     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 868     if (track->GetTPCNclsF()>0) {
 869       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 870     }
 871
 872     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 873
 874     eta=track->Eta();
 875     if ( TMath::Abs(eta)>0.9 ) continue;
 876
 877     //nacc++; // counter for accepted multiplicity
 878
 879     // the default ("basic") scenario
 880     if (scBasic == 1) {
 881       FillTPCHistogramsNsigma(fListQAtpcBasic,0,track,mult);
 882       FillTPCHistogramsSignal(fListQAtpcBasic,0,track,mult);
 883     }
 884
 885     // only MC truth identified particles
 886     if (scMCtruth == 1) {
 887       FillTPCHistogramsNsigma(fListQAtpcMCtruth,1,track,mult);
 888     }
 889
 890     // the "hybrid" scenario (only for LHC11h)
 891     if (scHybrid == 1) {
 892       FillTPCHistogramsNsigma(fListQAtpcHybrid,2,track,mult);
 893     }
 894
 895     // the "OROC high" scenario (only for LHC11h)
 896     if (scOROChigh == 1) {
 897       FillTPCHistogramsNsigma(fListQAtpcOROChigh,3,track,mult);
 898     }
 899
 900   } // -- End: track loop --
 901
 902
 903   // -- Begin: track loops for V0 candidates --
 904   if (scV0 == 1) {
 905
 906     // - Begin: track loop for electrons from V0 -
 907     for(Int_t itrack = 0; itrack < fV0electrons->GetEntries(); itrack++){
 908       AliVTrack *track=(AliVTrack*)fV0electrons->At(itrack);
 909
 910       //
 911       //basic track cuts
 912       //
 913       ULong_t status=track->GetStatus();
 914       // not that nice. status bits not in virtual interface
 915       // TPC refit + ITS refit + TPC pid
 916       if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 917           !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
 918
 919       // The TPC pid cut removes the light nuclei (>5 sigma from proton line)
 920       //||        !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  )
 921       Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 922       Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 923       if (track->GetTPCNclsF()>0) {
 924         ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 925       }
 926
 927       if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 928
 929       eta=track->Eta();
 930       if ( TMath::Abs(eta)>0.9 ) continue;
 931
 932       // fill histograms for V0 candidates
 933       FillTPCHistogramsNsigma(fListQAtpcV0,40,track,mult);
 934
 935     } // - End: track loop for electrons from V0 -
 936
 937
 938     // - Begin: track loop for pions from V0 -
 939     for(Int_t itrack = 0; itrack < fV0pions->GetEntries(); itrack++){
 940       AliVTrack *track=(AliVTrack*)fV0pions->At(itrack);
 941
 942       //
 943       //basic track cuts
 944       //
 945       ULong_t status=track->GetStatus();
 946       // not that nice. status bits not in virtual interface
 947       // TPC refit + ITS refit + TPC pid
 948       if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 949           !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
 950
 951       // The TPC pid cut removes the light nuclei (>5 sigma from proton line)
 952       //||        !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  )
 953       Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 954       Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 955       if (track->GetTPCNclsF()>0) {
 956         ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 957       }
 958
 959       if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 960
 961       eta=track->Eta();
 962       if ( TMath::Abs(eta)>0.9 ) continue;
 963
 964       // fill histograms for V0 candidates
 965       FillTPCHistogramsNsigma(fListQAtpcV0,42,track,mult);
 966
 967     } // - End: track loop for pions from V0 -
 968
 969
 970     // - Begin: track loop for kaons from V0 -
 971     for(Int_t itrack = 0; itrack < fV0kaons->GetEntries(); itrack++){
 972       AliVTrack *track=(AliVTrack*)fV0kaons->At(itrack);
 973
 974       //
 975       //basic track cuts
 976       //
 977       ULong_t status=track->GetStatus();
 978       // not that nice. status bits not in virtual interface
 979       // TPC refit + ITS refit + TPC pid
 980       if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 981           !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
 982
 983       // The TPC pid cut removes the light nuclei (>5 sigma from proton line)
 984       //||        !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  )
 985       Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 986       Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 987       if (track->GetTPCNclsF()>0) {
 988         ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 989       }
 990
 991       if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 992
 993       eta=track->Eta();
 994       if ( TMath::Abs(eta)>0.9 ) continue;
 995
 996       // fill histograms for V0 candidates
 997       FillTPCHistogramsNsigma(fListQAtpcV0,43,track,mult);
 998
 999     } // - End: track loop for kaons from V0 -
1000
1001
1002     // - Begin: track loop for protons from V0 -
1003     for(Int_t itrack = 0; itrack < fV0protons->GetEntries(); itrack++){
1004       AliVTrack *track=(AliVTrack*)fV0protons->At(itrack);
1005
1006       //
1007       //basic track cuts
1008       //
1009       ULong_t status=track->GetStatus();
1010       // not that nice. status bits not in virtual interface
1011       // TPC refit + ITS refit + TPC pid
1012       if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
1013           !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
1014
1015       // The TPC pid cut removes the light nuclei (>5 sigma from proton line)
1016       //||        !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  )
1017       Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
1018       Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1019       if (track->GetTPCNclsF()>0) {
1020         ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
1021       }
1022
1023       if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
1024
1025       eta=track->Eta();
1026       if ( TMath::Abs(eta)>0.9 ) continue;
1027
1028       // fill histograms for V0 candidates
1029       FillTPCHistogramsNsigma(fListQAtpcV0,44,track,mult);
1030
1031     } // - End: track loop for protons from V0 -
1032
1033   } // -- End: track loops for V0 candidates --
1034
1035
1036   // Multiplicity distribution
1037   TH1 *hm=(TH1*)fListQAtpc->At(scCounter);
1038   if (hm) {
1039     hm->Fill(mult);
1040   }
1041
1042   //printf("\nAccepted multiplicity = %i \n --- END of event --- \n",nacc);
1043
1044 }
1045
1046 //______________________________________________________________________________
1047 void AliAnalysisTaskPIDqa::FillTRDqa()
1048 {
1049   //
1050   // Fill PID qa histograms for the TRD
1051   //
1052   AliVEvent *event=InputEvent();
1053   Int_t ntracks = event->GetNumberOfTracks();
1054   for(Int_t itrack = 0; itrack <  ntracks; itrack++){
1055     AliVTrack *track = (AliVTrack *)event->GetTrack(itrack);
1056
1057     //
1058     //basic track cuts
1059     //
1060     ULong_t status=track->GetStatus();
1061     // not that nice. status bits not in virtual interface
1062     // TPC refit + ITS refit + TPC pid + TRD out
1063     if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
1064         !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
1065 //         !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  ) || //removes light nuclei. So it is out for the moment
1066         !( (status & AliVTrack::kTRDout  ) == AliVTrack::kTRDout  )) continue;
1067
1068     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
1069     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1070     if (track->GetTPCNclsF()>0) {
1071       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
1072     }
1073
1074     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
1075
1076     Double_t likelihoods[AliPID::kSPECIES];
1077     if(fPIDResponse->ComputeTRDProbability(track, AliPID::kSPECIES, likelihoods) != AliPIDResponse::kDetPidOk) continue;
1078     Int_t ntracklets = 0;
1079     Double_t momentum = -1.;
1080     for(Int_t itl = 0; itl < 6; itl++) {
1081       if(track->GetTRDmomentum(itl) > 0.) {
1082         ntracklets++;
1083         if(momentum < 0) momentum = track->GetTRDmomentum(itl);
1084       }
1085     }
1086
1087     for(Int_t ispecie = 0; ispecie < AliPID::kSPECIES; ispecie++){
1088       TH2F *hLike = (TH2F *)fListQAtrd->At(ntracklets*AliPID::kSPECIES+ispecie);
1089       if (hLike) hLike->Fill(momentum,likelihoods[ispecie]);
1090     }
1091
1092     //=== nSigma and signal ===
1093     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1094       TH2 *h=(TH2*)fListQAtrdNsig->At(ispecie);
1095       TH2 *hTPCTOF=(TH2*)fListQAtrdNsigTPCTOF->At(ispecie);
1096       if (!h || !hTPCTOF) continue;
1097       Float_t nSigmaTPC=fPIDResponse->NumberOfSigmas(AliPIDResponse::kTPC, track, (AliPID::EParticleType)ispecie);
1098       Float_t nSigmaTRD=fPIDResponse->NumberOfSigmas(AliPIDResponse::kTRD, track, (AliPID::EParticleType)ispecie);
1099       Float_t nSigmaTOF=fPIDResponse->NumberOfSigmas(AliPIDResponse::kTOF, track, (AliPID::EParticleType)ispecie);
1100       h->Fill(momentum,nSigmaTRD);
1101
1102       if (TMath::Abs(nSigmaTPC)<3 && TMath::Abs(nSigmaTOF)<3) {
1103         hTPCTOF->Fill(momentum,nSigmaTRD);
1104       }
1105     }
1106
1107     TH2 *h=(TH2*)fListQAtrdNsig->Last();
1108
1109     if (h) {
1110       Double_t sig=track->GetTRDsignal();
1111       h->Fill(momentum,sig);
1112     }
1113
1114   }
1115 }
1116
1117 //______________________________________________________________________________
1118 void AliAnalysisTaskPIDqa::FillTOFqa()
1119 {
1120   //
1121   // Fill TOF information
1122   //
1123   AliVEvent *event=InputEvent();
1124
1125   Int_t ntracks=event->GetNumberOfTracks();
1126   Int_t tracksAtTof = 0;
1127   for(Int_t itrack = 0; itrack < ntracks; itrack++){
1128     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
1129
1130     //
1131     //basic track cuts
1132     //
1133     ULong_t status=track->GetStatus();
1134     // TPC refit + ITS refit +
1135     // TOF out + kTIME
1136     // kTIME
1137     // (we don't use kTOFmismatch because it depends on TPC and kTOFpid because it prevents light nuclei
1138     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
1139         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
1140         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||
1141         //        !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||
1142         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;
1143
1144     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
1145     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1146     if (track->GetTPCNclsF()>0) {
1147       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
1148     }
1149
1150     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
1151
1152     tracksAtTof++;
1153
1154     Double_t mom=track->P();
1155
1156     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1157       TH2 *h=(TH2*)fListQAtof->At(ispecie);
1158       if (!h) continue;
1159       Double_t nSigma=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
1160       h->Fill(mom,nSigma);
1161     }
1162
1163     TH2 *h=(TH2*)fListQAtof->FindObject("hSigP_TOF");
1164     if (h) {
1165       Double_t sig=track->GetTOFsignal()/1000.;
1166       h->Fill(mom,sig);
1167     }
1168
1169     Int_t mask = fPIDResponse->GetTOFResponse().GetStartTimeMask(mom);
1170     ((TH1F*)fListQAtof->FindObject("hStartTimeMask_TOF"))->Fill((Double_t)(mask+0.5));
1171
1172     if (mom >= 0.75 && mom <= 1.25 ) {
1173       Double_t nsigma= fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)AliPID::kPion);
1174       if (mask == 0) {
1175         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-Fill"))->Fill(nsigma);
1176       } else if (mask == 1) {
1177         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-TOF"))->Fill(nsigma);
1178       } else if ( (mask == 2) || (mask == 4) || (mask == 6) ) {
1179         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-T0"))->Fill(nsigma);
1180       } else {
1181         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-Best"))->Fill(nsigma);
1182       }
1183       if (mask & 0x1) { //at least TOF-T0 present
1184         Double_t delta=0;
1185         (void)fPIDResponse->GetSignalDelta((AliPIDResponse::EDetector)AliPIDResponse::kTOF,track,(AliPID::EParticleType)AliPID::kPion,delta);
1186         ((TH1F*)fListQAtof->FindObject("hDelta_TOF_Pion"))->Fill(delta);
1187       }
1188     }
1189
1190     Double_t res = (Double_t)fPIDResponse->GetTOFResponse().GetStartTimeRes(mom);
1191     ((TH1F*)fListQAtof->FindObject("hStartTimeRes_TOF"))->Fill(res);
1192
1193     Double_t startTimeT0 = event->GetT0TOF(0);
1194     if (startTimeT0 < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeAC_T0"))->Fill(startTimeT0);
1195     else {
1196       startTimeT0 = event->GetT0TOF(1);
1197       if (startTimeT0 < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeA_T0"))->Fill(startTimeT0);
1198       startTimeT0 = event->GetT0TOF(2);
1199       if (startTimeT0 < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeC_T0"))->Fill(startTimeT0);
1200     }
1201   }
1202   if (tracksAtTof > 0) {
1203     ((TH1F* )fListQAtof->FindObject("hnTracksAt_TOF"))->Fill(tracksAtTof);
1204     Int_t mask = fPIDResponse->GetTOFResponse().GetStartTimeMask(5.);
1205     if (mask & 0x1) ((TH1F*)fListQAtof->FindObject("hT0MakerEff"))->Fill(tracksAtTof);
1206   }
1207 }
1208
1209 //______________________________________________________________________________
1210 void AliAnalysisTaskPIDqa::FillT0qa()
1211 {
1212   //
1213   // Fill TOF information
1214   //
1215   AliVEvent *event=InputEvent();
1216
1217   Int_t ntracks=event->GetNumberOfTracks();
1218
1219   Int_t tracksAtT0 = 0;
1220
1221   for(Int_t itrack = 0; itrack < ntracks; itrack++){
1222     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
1223
1224     //
1225     //basic track cuts
1226     //
1227     ULong_t status=track->GetStatus();
1228     // TPC refit + ITS refit +
1229     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
1230         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
1231     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
1232     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1233     if (track->GetTPCNclsF()>0) {
1234       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
1235     }
1236     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
1237
1238     tracksAtT0++;
1239   }
1240
1241   Bool_t t0A = kFALSE;
1242   Bool_t t0C = kFALSE;
1243   Bool_t t0And = kFALSE;
1244   Double_t startTimeT0 = event->GetT0TOF(0);     // AND
1245   if (startTimeT0 < 90000) {
1246     t0And = kTRUE;
1247     ((TH1F*)fListQAt0->FindObject("hStartTimeAC_T0"))->Fill(startTimeT0);
1248     }
1249   startTimeT0 = event->GetT0TOF(1);             // T0A
1250   if (startTimeT0 < 90000) {
1251     t0A = kTRUE;
1252     ((TH1F*)fListQAt0->FindObject("hStartTimeA_T0"))->Fill(startTimeT0);
1253
1254   }
1255   startTimeT0 = event->GetT0TOF(2);             // T0C
1256   if (startTimeT0 < 90000) {
1257     t0C = kTRUE;
1258     ((TH1F*)fListQAt0->FindObject("hStartTimeC_T0"))->Fill(startTimeT0);
1259   }
1260
1261   ((TH1F* )fListQAt0->FindObject("hnTracksAt_T0"))->Fill(tracksAtT0);
1262   if (t0A) ((TH1F*)fListQAt0->FindObject("hT0AEff"))->Fill(tracksAtT0);
1263   if (t0C) ((TH1F*)fListQAt0->FindObject("hT0CEff"))->Fill(tracksAtT0);
1264   if (t0And) ((TH1F*)fListQAt0->FindObject("hT0AndEff"))->Fill(tracksAtT0);
1265   if (t0A || t0C) ((TH1F*)fListQAt0->FindObject("hT0OrEff"))->Fill(tracksAtT0);
1266 }
1267
1268
1269 //______________________________________________________________________________
1270 void AliAnalysisTaskPIDqa::FillEMCALqa()
1271 {
1272   //
1273   // Fill PID qa histograms for the EMCAL
1274   //
1275
1276   AliVEvent *event=InputEvent();
1277
1278   Int_t ntracks=event->GetNumberOfTracks();
1279   for(Int_t itrack = 0; itrack < ntracks; itrack++){
1280     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
1281
1282     //
1283     //basic track cuts
1284     //
1285     ULong_t status=track->GetStatus();
1286     // not that nice. status bits not in virtual interface
1287     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
1288
1289     Double_t pt=track->Pt();
1290
1291     //EMCAL nSigma (only for electrons at the moment)
1292     TH2 *h=(TH2*)fListQAemcal->At(0);
1293     if (!h) continue;
1294     Double_t nSigma=fPIDResponse->NumberOfSigmasEMCAL(track, (AliPID::EParticleType)0);
1295     h->Fill(pt,nSigma);
1296
1297   }
1298
1299    //EMCAL signal (E/p vs. pT) for electrons from V0
1300   for(Int_t itrack = 0; itrack < fV0electrons->GetEntries(); itrack++){
1301     AliVTrack *track=(AliVTrack*)fV0electrons->At(itrack);
1302
1303     //
1304     //basic track cuts
1305     //
1306     ULong_t status=track->GetStatus();
1307     // not that nice. status bits not in virtual interface
1308     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
1309
1310     Double_t pt=track->Pt();
1311
1312     TH2 *h=(TH2*)fListQAemcal->At(1);
1313     if (h) {
1314
1315       Int_t nMatchClus = track->GetEMCALcluster();
1316       Double_t mom     = track->P();
1317       Double_t eop     = -1.;
1318
1319       if(nMatchClus > -1){
1320
1321         AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
1322
1323         if(matchedClus){
1324
1325           // matched cluster is EMCAL
1326           if(matchedClus->IsEMCAL()){
1327
1328             Double_t fClsE       = matchedClus->E();
1329             eop                  = fClsE/mom;
1330
1331             h->Fill(pt,eop);
1332
1333           }
1334         }
1335       }
1336     }
1337   }
1338
1339    //EMCAL signal (E/p vs. pT) for pions from V0
1340   for(Int_t itrack = 0; itrack < fV0pions->GetEntries(); itrack++){
1341     AliVTrack *track=(AliVTrack*)fV0pions->At(itrack);
1342
1343     //
1344     //basic track cuts
1345     //
1346     ULong_t status=track->GetStatus();
1347     // not that nice. status bits not in virtual interface
1348     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
1349
1350     Double_t pt=track->Pt();
1351
1352     TH2 *h=(TH2*)fListQAemcal->At(2);
1353     if (h) {
1354
1355       Int_t nMatchClus = track->GetEMCALcluster();
1356       Double_t mom     = track->P();
1357       Double_t eop     = -1.;
1358
1359       if(nMatchClus > -1){
1360
1361         AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
1362
1363         if(matchedClus){
1364
1365           // matched cluster is EMCAL
1366           if(matchedClus->IsEMCAL()){
1367
1368             Double_t fClsE       = matchedClus->E();
1369             eop                  = fClsE/mom;
1370
1371             h->Fill(pt,eop);
1372
1373           }
1374         }
1375       }
1376     }
1377   }
1378
1379    //EMCAL signal (E/p vs. pT) for protons from V0
1380   for(Int_t itrack = 0; itrack < fV0protons->GetEntries(); itrack++){
1381     AliVTrack *track=(AliVTrack*)fV0protons->At(itrack);
1382
1383     //
1384     //basic track cuts
1385     //
1386     ULong_t status=track->GetStatus();
1387     // not that nice. status bits not in virtual interface
1388     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
1389
1390     Double_t pt=track->Pt();
1391
1392     TH2 *hP=(TH2*)fListQAemcal->At(3);
1393     TH2 *hAP=(TH2*)fListQAemcal->At(4);
1394     if (hP && hAP) {
1395
1396       Int_t nMatchClus = track->GetEMCALcluster();
1397       Double_t mom     = track->P();
1398       Int_t charge     = track->Charge();
1399       Double_t eop     = -1.;
1400
1401       if(nMatchClus > -1){
1402
1403         AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
1404
1405         if(matchedClus){
1406
1407           // matched cluster is EMCAL
1408           if(matchedClus->IsEMCAL()){
1409
1410             Double_t fClsE       = matchedClus->E();
1411             eop                  = fClsE/mom;
1412
1413             if(charge > 0)      hP->Fill(pt,eop);
1414             else if(charge < 0) hAP->Fill(pt,eop);
1415
1416           }
1417         }
1418       }
1419     }
1420   }
1421
1422 }
1423
1424
1425 //______________________________________________________________________________
1426 void AliAnalysisTaskPIDqa::FillHMPIDqa()
1427 {
1428   //
1429   // Fill PID qa histograms for the HMPID
1430   //
1431
1432   AliVEvent *event=InputEvent();
1433
1434   Int_t ntracks=event->GetNumberOfTracks();
1435   for(Int_t itrack = 0; itrack < ntracks; itrack++){
1436     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
1437
1438     //
1439     //basic track cuts
1440     //
1441     const ULong_t status=track->GetStatus();
1442     // not that nice. status bits not in virtual interface
1443     // TPC refit + ITS refit +
1444     // TOF out + TOFpid +
1445     // kTIME
1446     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
1447         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
1448
1449     const Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
1450     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1451     if (track->GetTPCNclsF()>0) {
1452       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
1453     }
1454
1455     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
1456
1457     const Double_t mom = track->P();
1458     const Double_t ckovAngle = track->GetHMPIDsignal();
1459
1460     Int_t nhists=0;
1461     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1462       if (ispecie==AliPID::kElectron || ispecie==AliPID::kMuon) continue;
1463       TH2 *h=(TH2*)fListQAhmpid->At(nhists);
1464       if (!h) {++nhists; continue;}
1465       const Double_t nSigma=fPIDResponse->NumberOfSigmasHMPID(track, (AliPID::EParticleType)ispecie);
1466       h->Fill(mom,nSigma);
1467       ++nhists;
1468     }
1469
1470     TH1F *hThetavsMom = (TH1F*)fListQAhmpid->At(AliPID::kSPECIESC);
1471
1472     if (hThetavsMom) hThetavsMom->Fill(mom,ckovAngle);
1473
1474   }
1475 }
1476 //______________________________________________________________________________
1477 void AliAnalysisTaskPIDqa::FillTOFHMPIDqa()
1478 {
1479   //
1480   // Fill PID qa histograms for the HMPID
1481   //
1482
1483   AliVEvent *event=InputEvent();
1484
1485   Int_t ntracks=event->GetNumberOfTracks();
1486   for(Int_t itrack = 0; itrack < ntracks; itrack++){
1487     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
1488
1489     //
1490     //basic track cuts
1491     //
1492     ULong_t status=track->GetStatus();
1493     // not that nice. status bits not in virtual interface
1494     // TPC refit + ITS refit +
1495     // TOF out + TOFpid +
1496     // kTIME
1497     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
1498         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
1499         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||
1500         !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||
1501         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;
1502
1503     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
1504     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1505     if (track->GetTPCNclsF()>0) {
1506       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
1507     }
1508
1509     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
1510
1511     Double_t mom = track->P();
1512     Double_t ckovAngle = track->GetHMPIDsignal();
1513
1514     Double_t nSigmaTOF[3];
1515     TH1F *h[3];
1516
1517     for (Int_t ispecie=2; ispecie<5; ++ispecie){
1518       //TOF nSigma
1519       nSigmaTOF[ispecie-2]=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
1520       h[ispecie-2] = (TH1F*)fListQAtofhmpid->At(ispecie-2);}
1521
1522     if(TMath::Abs(nSigmaTOF[0])<2)                                                              h[0]->Fill(mom,ckovAngle);
1523
1524     if(TMath::Abs(nSigmaTOF[1])<2 && TMath::Abs(nSigmaTOF[0])>3)                                h[1]->Fill(mom,ckovAngle);
1525
1526     if(TMath::Abs(nSigmaTOF[2])<2 && TMath::Abs(nSigmaTOF[1])>3 && TMath::Abs(nSigmaTOF[0])>3)  h[2]->Fill(mom,ckovAngle);
1527
1528   }
1529
1530 }
1531
1532 //______________________________________________________________________________
1533 void AliAnalysisTaskPIDqa::FillTPCTOFqa()
1534 {
1535   //
1536   // Fill PID qa histograms for the TOF
1537   //   Here also the TPC histograms after TOF selection are filled
1538   //
1539
1540   AliVEvent *event=InputEvent();
1541
1542   Int_t ntracks=event->GetNumberOfTracks();
1543   for(Int_t itrack = 0; itrack < ntracks; itrack++){
1544     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
1545
1546     //
1547     //basic track cuts
1548     //
1549     ULong_t status=track->GetStatus();
1550     // not that nice. status bits not in virtual interface
1551     // TPC refit + ITS refit +
1552     // TOF out + TOFpid +
1553     // kTIME
1554     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
1555         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
1556 //         !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid ) || //removes light nuclei, so it is out for the moment
1557         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||
1558         !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||
1559         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;
1560
1561     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
1562     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
1563     if (track->GetTPCNclsF()>0) {
1564       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
1565     }
1566
1567     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
1568
1569
1570     Double_t mom=track->P();
1571     Double_t momTPC=track->GetTPCmomentum();
1572
1573     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1574       //TOF nSigma
1575       Double_t nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
1576       Double_t nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
1577
1578       //TPC after TOF cut
1579       TH2 *h=(TH2*)fListQAtpctof->At(ispecie);
1580       if (h && TMath::Abs(nSigmaTOF)<3.) h->Fill(momTPC,nSigmaTPC);
1581
1582       //TOF after TPC cut
1583       h=(TH2*)fListQAtpctof->At(ispecie+AliPID::kSPECIESC);
1584       if (h && TMath::Abs(nSigmaTPC)<3.) h->Fill(mom,nSigmaTOF);
1585
1586       //EMCAL after TOF and TPC cut
1587       h=(TH2*)fListQAtpctof->At(ispecie+2*AliPID::kSPECIESC);
1588       if (h && TMath::Abs(nSigmaTOF)<3. && TMath::Abs(nSigmaTPC)<3. ){
1589
1590         Int_t nMatchClus = track->GetEMCALcluster();
1591         Double_t pt      = track->Pt();
1592         Double_t eop     = -1.;
1593
1594         if(nMatchClus > -1){
1595
1596           AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
1597
1598           if(matchedClus){
1599
1600             // matched cluster is EMCAL
1601             if(matchedClus->IsEMCAL()){
1602
1603               Double_t fClsE       = matchedClus->E();
1604               eop                  = fClsE/mom;
1605
1606               h->Fill(pt,eop);
1607
1608
1609             }
1610           }
1611         }
1612       }
1613     }
1614   }
1615 }
1616
1617 //_____________________________________________________________________________
1618 void AliAnalysisTaskPIDqa::FillQAinfo()
1619 {
1620   //
1621   // Fill the QA information
1622   //
1623
1624
1625   //TPC QA info
1626   TObjArray *arrTPC=static_cast<TObjArray*>(fListQAinfo->At(0));
1627   if (fPIDResponse && arrTPC){
1628     AliTPCPIDResponse &tpcResp=fPIDResponse->GetTPCResponse();
1629     // fill spline names
1630     if (!arrTPC->UncheckedAt(0)){
1631
1632       TObjArray *arrTPCsplineNames=new TObjArray(AliPID::kSPECIESC);
1633       arrTPCsplineNames->SetOwner();
1634       arrTPCsplineNames->SetName("TPC_spline_names");
1635       arrTPC->AddAt(arrTPCsplineNames,0);
1636
1637       for (Int_t iresp=0; iresp<AliPID::kSPECIESC; ++iresp){
1638         const TObject *o=tpcResp.GetResponseFunction((AliPID::EParticleType)iresp);
1639         if (!o) continue;
1640         arrTPCsplineNames->Add(new TObjString(Form("%02d: %s",iresp, o->GetName())));
1641       }
1642     }
1643
1644     // tpc response config
1645     if (!arrTPC->UncheckedAt(1)){
1646
1647       TObjArray *arrTPCconfigInfo=new TObjArray;
1648       arrTPCconfigInfo->SetOwner();
1649       arrTPCconfigInfo->SetName("TPC_config_info");
1650       arrTPC->AddAt(arrTPCconfigInfo,1);
1651
1652       TObjString *ostr=0x0;
1653       ostr=new TObjString;
1654       ostr->String().Form("Eta Corr map: %s", tpcResp.GetEtaCorrMap()?tpcResp.GetEtaCorrMap()->GetName():"none");
1655       arrTPCconfigInfo->Add(ostr);
1656
1657       ostr=new TObjString;
1658       ostr->String().Form("Sigma Par map: %s", tpcResp.GetSigmaPar1Map()?tpcResp.GetSigmaPar1Map()->GetName():"none");
1659       arrTPCconfigInfo->Add(ostr);
1660
1661       ostr=new TObjString;
1662       ostr->String().Form("MIP: %.2f", tpcResp.GetMIP());
1663       arrTPCconfigInfo->Add(ostr);
1664
1665       ostr=new TObjString;
1666       ostr->String().Form("Res: Def %.3g (%.3g) : AllHigh %.3g (%.3g) : OROC high %.3g (%.3g)",
1667                           tpcResp.GetRes0(AliTPCPIDResponse::kDefault), tpcResp.GetResN2(AliTPCPIDResponse::kDefault),
1668                           tpcResp.GetRes0(AliTPCPIDResponse::kALLhigh), tpcResp.GetResN2(AliTPCPIDResponse::kALLhigh),
1669                           tpcResp.GetRes0(AliTPCPIDResponse::kOROChigh), tpcResp.GetResN2(AliTPCPIDResponse::kOROChigh)
1670                          );
1671       arrTPCconfigInfo->Add(ostr);
1672     }
1673   }
1674 }
1675
1676 //______________________________________________________________________________
1677 void AliAnalysisTaskPIDqa::SetupITSqa()
1678 {
1679   //
1680   // Create the ITS qa objects
1681   //
1682
1683   TVectorD *vX=MakeLogBinning(200,.1,30);
1684
1685   //ITS+TPC tracks
1686   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1687     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_ITS_%s",AliPID::ParticleName(ispecie)),
1688                               Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1689                               vX->GetNrows()-1,vX->GetMatrixArray(),
1690                               200,-10,10);
1691     fListQAits->Add(hNsigmaP);
1692   }
1693   TH2F *hSig = new TH2F("hSigP_ITS",
1694                         "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
1695                         vX->GetNrows()-1,vX->GetMatrixArray(),
1696                         300,0,300);
1697   fListQAits->Add(hSig);
1698
1699   //ITS Standalone tracks
1700   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1701     TH2F *hNsigmaPSA = new TH2F(Form("hNsigmaP_ITSSA_%s",AliPID::ParticleName(ispecie)),
1702                                 Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1703                                 vX->GetNrows()-1,vX->GetMatrixArray(),
1704                                 200,-10,10);
1705     fListQAitsSA->Add(hNsigmaPSA);
1706   }
1707   TH2F *hSigSA = new TH2F("hSigP_ITSSA",
1708                           "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
1709                           vX->GetNrows()-1,vX->GetMatrixArray(),
1710                           300,0,300);
1711   fListQAitsSA->Add(hSigSA);
1712
1713   //ITS Pure Standalone tracks
1714   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1715     TH2F *hNsigmaPPureSA = new TH2F(Form("hNsigmaP_ITSPureSA_%s",AliPID::ParticleName(ispecie)),
1716                                     Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1717                                     vX->GetNrows()-1,vX->GetMatrixArray(),
1718                                     200,-10,10);
1719     fListQAitsPureSA->Add(hNsigmaPPureSA);
1720   }
1721   TH2F *hSigPureSA = new TH2F("hSigP_ITSPureSA",
1722                               "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
1723                               vX->GetNrows()-1,vX->GetMatrixArray(),
1724                               300,0,300);
1725   fListQAitsPureSA->Add(hSigPureSA);
1726
1727   delete vX;
1728 }
1729
1730 //_____________________________________________________________________________
1731 void AliAnalysisTaskPIDqa::AddTPCHistogramsSignal(TList *sublist, const char *scenario)
1732 {
1733   //
1734   // Create the TPC qa objects: create histograms for the TPC signal for different settings
1735   //
1736
1737   TVectorD *vX=MakeLogBinning(200,.1,30);
1738   Int_t nBinsMult = 38;
1739   Double_t xBinsMult[39] = {0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100,
1740                            120, 140, 160, 180, 200,
1741                            300, 400, 500, 600, 700, 800, 900, 1000,
1742                            1200, 1400, 1600, 1800, 2000,
1743                            2200, 2400, 2600, 2800, 3000,
1744                            3200, 3400, 3600, 3800, 4000
1745                            };
1746   const Int_t binsEta=110;
1747   Float_t etaMin=-1.1;
1748   Float_t etaMax=1.1;
1749
1750   char signal[4][12]={"std","IROC","OROCmedium","OROClong"};
1751
1752
1753   // TPC signal vs. p for all particles (standard, IROC, OROCmedium, OROClong)
1754   for (Int_t iSig=0; iSig<4; iSig++) {
1755     TH2F *hSigP = new TH2F(Form("hSigP_TPC_%s_%s",signal[iSig],scenario),
1756                              Form("TPC_%s signal (%s) vs. p;p [GeV]; TPC signal [arb. units]",scenario,signal[iSig]),
1757                              vX->GetNrows()-1,vX->GetMatrixArray(),
1758                              300,0,300);
1759     sublist->Add(hSigP);
1760   }
1761
1762   // MIP pions: TPC signal vs. eta
1763   for (Int_t iSig=0; iSig<4; iSig++) {
1764     TH2F *hSigEtaMIPpi = new TH2F(Form("hSigEta_TPC_%s_%s_MIPpi",signal[iSig],scenario),
1765                                   Form("TPC_%s signal (%s) MIPpi vs. eta;#eta;TPC signal [arb. units]",scenario,signal[iSig]),
1766                                   binsEta,etaMin,etaMax,
1767                                   300,0,300);
1768     sublist->Add(hSigEtaMIPpi);
1769   }
1770
1771   // MIP pions: TPC signal vs. multiplicity
1772   for (Int_t iSig=0; iSig<4; iSig++) {
1773     TH2F *hSigMultMPIpi = new TH2F(Form("hSigMult_TPC_%s_%s_MIPpi",signal[iSig],scenario),
1774                                    Form("TPC_%s signal (%s) MIPpi vs. mult;multiplicity;TPC signal [arb. units]",scenario,signal[iSig]),
1775                                    nBinsMult,xBinsMult,
1776                                    300,0,300);
1777     sublist->Add(hSigMultMPIpi);
1778   }
1779
1780   // Electrons: TPC signal vs. eta
1781   for (Int_t iSig=0; iSig<4; iSig++) {
1782     TH2F *hSigEtaEle = new TH2F(Form("hSigEta_TPC_%s_%s_Ele",signal[iSig],scenario),
1783                                 Form("TPC_%s signal (%s) electrons vs. eta;#eta;TPC signal [arb. units]",scenario,signal[iSig]),
1784                                 binsEta,etaMin,etaMax,
1785                                 300,0,300);
1786     sublist->Add(hSigEtaEle);
1787   }
1788
1789   // Electrons: TPC signal vs. multiplicity
1790   for (Int_t iSig=0; iSig<4; iSig++) {
1791     TH2F *hSigMultEle = new TH2F(Form("hSigMult_TPC_%s_%s_Ele",signal[iSig],scenario),
1792                                  Form("TPC_%s signal (%s) electrons vs. mult;multiplicity;TPC signal [arb. units]",scenario,signal[iSig]),
1793                                  nBinsMult,xBinsMult,
1794                                  300,0,300);
1795     sublist->Add(hSigMultEle);
1796   }
1797
1798   delete vX;
1799
1800 }
1801
1802 //_____________________________________________________________________________
1803 void AliAnalysisTaskPIDqa::AddTPCHistogramsNsigma(TList *sublist, const char *scenario, Int_t scnumber)
1804 {
1805   //
1806   // Create the TPC qa objects: create histograms for TPC Nsigma for different settings
1807   //
1808
1809   TVectorD *vX=MakeLogBinning(200,.1,30.);
1810   Int_t nBinsMult = 38;
1811   Double_t xBinsMult[39] = {0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100,
1812                            120, 140, 160, 180, 200,
1813                            300, 400, 500, 600, 700, 800, 900, 1000,
1814                            1200, 1400, 1600, 1800, 2000,
1815                            2200, 2400, 2600, 2800, 3000,
1816                            3200, 3400, 3600, 3800, 4000
1817                            };
1818   const Int_t binsEta=110;
1819   Float_t etaMin=-1.1;
1820   Float_t etaMax=1.1;
1821
1822   Int_t nSpecies=0;
1823
1824   if (scnumber == 4) nSpecies=(Int_t)AliPID::kSPECIES;
1825   else nSpecies=(Int_t)AliPID::kSPECIESC;
1826
1827   // Nsigma vs. p for different particle species
1828   for (Int_t ispecie=0; ispecie<nSpecies; ++ispecie){
1829     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_%s_%s",scenario,AliPID::ParticleName(ispecie)),
1830                               Form("TPC_%s n#sigma %s vs. p;p [GeV]; n#sigma",scenario,AliPID::ParticleName(ispecie)),
1831                               vX->GetNrows()-1,vX->GetMatrixArray(),
1832                               200,-10,10);
1833     sublist->Add(hNsigmaP);
1834   }
1835
1836   // Nsigma vs. eta for different particle species (only for some scenarios)
1837   if ( scnumber == 1 || scnumber == 4 ) {
1838     for (Int_t ispecie=0; ispecie<nSpecies; ++ispecie){
1839       TH2F *hNsigmaEta = new TH2F(Form("hNsigmaEta_TPC_%s_%s",scenario,AliPID::ParticleName(ispecie)),
1840                               Form("TPC_%s n#sigma %s vs. eta;#eta; n#sigma",scenario,AliPID::ParticleName(ispecie)),
1841                               binsEta,etaMin,etaMax,
1842                               200,-10,10);
1843       sublist->Add(hNsigmaEta);
1844     }
1845   }
1846
1847   // Nsigma vs. multiplicity for different particle species (only for some scenarios)
1848   if ( scnumber == 1 || scnumber == 4 ) {
1849     for (Int_t ispecie=0; ispecie<nSpecies; ++ispecie){
1850       TH2F *hNsigmaMult = new TH2F(Form("hNsigmaMult_TPC_%s_%s",scenario,AliPID::ParticleName(ispecie)),
1851                               Form("TPC_%s n#sigma %s vs. mult;multiplicity; n#sigma",scenario,AliPID::ParticleName(ispecie)),
1852                               nBinsMult,xBinsMult,
1853                               200,-10,10);
1854       sublist->Add(hNsigmaMult);
1855     }
1856   }
1857
1858   // - Beginn: Adding histograms for MIP pions and electrons (only for some scenarios) -
1859   if ( scnumber == 0 || scnumber == 2 || scnumber == 3 ) {
1860
1861     // MIP pions: Nsigma vs. eta
1862     TH2F *hNsigmaEtaMIPpi = new TH2F(Form("hNsigmaEta_TPC_%s_MIPpi",scenario),
1863                               Form("TPC_%s n#sigma MIPpi vs. eta;#eta; n#sigma",scenario),
1864                               binsEta,etaMin,etaMax,
1865                               200,-10,10);
1866     sublist->Add(hNsigmaEtaMIPpi);
1867
1868     // MIP pions: Nsigma vs. multiplicity
1869     TH2F *hNsigmaMultMIPpi = new TH2F(Form("hNsigmaMult_TPC_%s_MIPpi",scenario),
1870                                Form("TPC_%s n#sigma MIPpi vs. mult;multiplicity; n#sigma",scenario),
1871                                nBinsMult,xBinsMult,
1872                                200,-10,10);
1873     sublist->Add(hNsigmaMultMIPpi);
1874
1875     // Electrons: Nsigma vs. eta
1876     TH2F *hNsigmaEtaEle = new TH2F(Form("hNsigmaEta_TPC_%s_Ele",scenario),
1877                               Form("TPC_%s n#sigma electrons vs. eta;#eta; n#sigma",scenario),
1878                               binsEta,etaMin,etaMax,
1879                               200,-10,10);
1880     sublist->Add(hNsigmaEtaEle);
1881
1882     // Electrons: Nsigma vs. multiplicity
1883     TH2F *hNsigmaMultEle = new TH2F(Form("hNsigmaMult_TPC_%s_Ele",scenario),
1884                                Form("TPC_%s n#sigma electrons vs. mult;multiplicity; n#sigma",scenario),
1885                                nBinsMult,xBinsMult,
1886                                200,-10,10);
1887     sublist->Add(hNsigmaMultEle);
1888   } // - End: Adding histograms for MIP pions and electrons
1889
1890   delete vX;
1891
1892 }
1893
1894 //______________________________________________________________________________
1895 void AliAnalysisTaskPIDqa::SetupTPCqa(Bool_t fillMC, Bool_t fill11h, Bool_t fillV0)
1896 {
1897   //
1898   // Create the TPC qa objects
1899   //
1900
1901   // Set up the multiplicity binning
1902   Int_t nBinsMult = 38;
1903   Double_t xBinsMult[39] = {0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100,
1904                            120, 140, 160, 180, 200,
1905                            300, 400, 500, 600, 700, 800, 900, 1000,
1906                            1200, 1400, 1600, 1800, 2000,
1907                            2200, 2400, 2600, 2800, 3000,
1908                            3200, 3400, 3600, 3800, 4000
1909                            };
1910
1911
1912   // Create TPC sublists for different scenarios
1913   // corresponding to available information,
1914   // e.g. MC or not, special settings for LHC11h
1915
1916   // basic/default scenario, used always
1917   fListQAtpcBasic=new TList;
1918   fListQAtpcBasic->SetOwner();
1919   fListQAtpcBasic->SetName("TPCBasic");
1920   fListQAtpc->Add(fListQAtpcBasic);
1921
1922   // MC truth scenario: use only MC truth identified particles
1923   // only available for MC
1924   if (fillMC == kTRUE) {
1925     fListQAtpcMCtruth=new TList;
1926     fListQAtpcMCtruth->SetOwner();
1927     fListQAtpcMCtruth->SetName("TPCMCtruth");
1928     fListQAtpc->Add(fListQAtpcMCtruth);
1929   }
1930
1931   // Hybrid and OROChigh scenarios,
1932   // special settings only available for PbPb LHC11h data
1933   if (fill11h == kTRUE) {
1934     fListQAtpcHybrid=new TList;
1935     fListQAtpcHybrid->SetOwner();
1936     fListQAtpcHybrid->SetName("TPCHybrid");
1937     fListQAtpc->Add(fListQAtpcHybrid);
1938
1939     fListQAtpcOROChigh=new TList;
1940     fListQAtpcOROChigh->SetOwner();
1941     fListQAtpcOROChigh->SetName("TPCOROChigh");
1942     fListQAtpc->Add(fListQAtpcOROChigh);
1943   }
1944
1945   // scenario only for V0s,
1946   // only available for ESDs
1947   if (fillV0 == kTRUE) {
1948     fListQAtpcV0=new TList;
1949     fListQAtpcV0->SetOwner();
1950     fListQAtpcV0->SetName("TPCV0");
1951     fListQAtpc->Add(fListQAtpcV0);
1952   }
1953
1954
1955   // the default ("basic") scenario
1956   AddTPCHistogramsNsigma(fListQAtpcBasic,"Basic",0);
1957   AddTPCHistogramsSignal(fListQAtpcBasic,"Basic");
1958
1959   // only MC truth identified particles
1960   if (fillMC) {
1961     AddTPCHistogramsNsigma(fListQAtpcMCtruth,"MCtruth",1);
1962   }
1963
1964   // the "hybrid" scenario (only for period LHC11h)
1965   if (fill11h) {
1966     AddTPCHistogramsNsigma(fListQAtpcHybrid,"Hybrid",2);
1967   }
1968
1969   // the "OROC high" scenario (only for period LHC11h)
1970   if (fill11h) {
1971     AddTPCHistogramsNsigma(fListQAtpcOROChigh,"OROChigh",3);
1972   }
1973
1974   // only for V0s
1975   if (fillV0) {
1976     AddTPCHistogramsNsigma(fListQAtpcV0,"V0",4);
1977   }
1978
1979
1980   // Multiplicity distribution --- as check
1981   TH1F *hMult = new TH1F("hMult_TPC",
1982                          "Multiplicity distribution;multiplicity;counts",
1983                          nBinsMult,xBinsMult);
1984   fListQAtpc->Add(hMult);
1985
1986 }
1987
1988 //______________________________________________________________________________
1989 void AliAnalysisTaskPIDqa::SetupTRDqa()
1990 {
1991   //
1992   // Create the TRD qa objects
1993   //
1994   TVectorD *vX=MakeLogBinning(200,.1,30);
1995   for(Int_t itl = 0; itl < 6; ++itl){
1996     for(Int_t ispecie = 0; ispecie < AliPID::kSPECIES; ispecie++){
1997       TH2F *hLikeP = new TH2F(Form("hLikeP_TRD_%dtls_%s", itl, AliPID::ParticleName(ispecie)),
1998                               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)),
1999                               vX->GetNrows()-1, vX->GetMatrixArray(),
2000                               100, 0., 1.);
2001       fListQAtrd->Add(hLikeP);
2002     }
2003   }
2004
2005   // === nSigma Values and signal ===
2006   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
2007     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TRD_%s",AliPID::ParticleName(ispecie)),
2008                               Form("TRD n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
2009                               vX->GetNrows()-1,vX->GetMatrixArray(),
2010                               100,-10,10);
2011     fListQAtrdNsig->Add(hNsigmaP);
2012   }
2013
2014   TH2F *hSig = new TH2F("hSigP_TRD",
2015                         "TRD signal vs. p;p [GeV]; TRD signal [arb. units]",
2016                         vX->GetNrows()-1,vX->GetMatrixArray(),
2017                         100,0,100);
2018   fListQAtrdNsig->Add(hSig);
2019
2020   fListQAtrd->Add(fListQAtrdNsig);
2021
2022   // === Same after 3 sigma in TPC and TOF
2023   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
2024     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TRD_TPCTOF_%s",AliPID::ParticleName(ispecie)),
2025                               Form("TRD n#sigma %s vs. p after 3#sigma cut in TPC&TOF;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
2026                               vX->GetNrows()-1,vX->GetMatrixArray(),
2027                               100,-10,10);
2028     fListQAtrdNsigTPCTOF->Add(hNsigmaP);
2029   }
2030
2031   fListQAtrd->Add(fListQAtrdNsigTPCTOF);
2032
2033   delete vX;
2034 }
2035
2036 //______________________________________________________________________________
2037 void AliAnalysisTaskPIDqa::SetupTOFqa()
2038 {
2039   //
2040   // Create the TOF qa objects
2041   //
2042
2043   TVectorD *vX=MakeLogBinning(200,.1,30);
2044
2045   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
2046     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_%s",AliPID::ParticleName(ispecie)),
2047                               Form("TOF n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
2048                               vX->GetNrows()-1,vX->GetMatrixArray(),
2049                               200,-10,10);
2050     fListQAtof->Add(hNsigmaP);
2051   }
2052
2053   TH1F *hnSigT0Fill = new TH1F("hNsigma_TOF_Pion_T0-Fill","TOF n#sigma (Pion) T0-FILL [0.75-1.25. GeV/c]",200,-10,10);
2054   fListQAtof->Add(hnSigT0Fill);
2055   TH1F *hnSigT0T0 = new TH1F("hNsigma_TOF_Pion_T0-T0","TOF n#sigma (Pion) T0-T0 [0.75-1.25 GeV/c]",200,-10,10);
2056   fListQAtof->Add(hnSigT0T0);
2057   TH1F *hnSigT0TOF = new TH1F("hNsigma_TOF_Pion_T0-TOF","TOF n#sigma (Pion) T0-TOF [0.75-1.25 GeV/c]",200,-10,10);
2058   fListQAtof->Add(hnSigT0TOF);
2059   TH1F *hnSigT0Best = new TH1F("hNsigma_TOF_Pion_T0-Best","TOF n#sigma (Pion) T0-Best [0.75-1.25 GeV/c]",200,-10,10);
2060   fListQAtof->Add(hnSigT0Best);
2061   TH1F *hnDeltaPi = new TH1F("hDelta_TOF_Pion","DeltaT (Pion) [0.75-1.25 GeV/c]",50,-500,500);
2062   fListQAtof->Add(hnDeltaPi);
2063
2064   TH2F *hSig = new TH2F("hSigP_TOF",
2065                         "TOF signal vs. p;p [GeV]; TOF signal [ns]",
2066                         vX->GetNrows()-1,vX->GetMatrixArray(),
2067                         300,0,30);
2068
2069   delete vX;
2070
2071   fListQAtof->Add(hSig);
2072
2073   TH1F *hStartTimeMaskTOF = new TH1F("hStartTimeMask_TOF","StartTime mask",8,0,8);
2074   fListQAtof->Add(hStartTimeMaskTOF);
2075   TH1F *hStartTimeResTOF = new TH1F("hStartTimeRes_TOF","StartTime resolution [ps]",100,0,500);
2076   fListQAtof->Add(hStartTimeResTOF);
2077
2078   TH1F *hnTracksAtTOF = new TH1F("hnTracksAt_TOF","Matched tracks at TOF",100,0,100);
2079   fListQAtof->Add(hnTracksAtTOF);
2080   TH1F *hT0MakerEff = new TH1F("hT0MakerEff","Events with T0-TOF vs nTracks",100,0,100);
2081   fListQAtof->Add(hT0MakerEff);
2082
2083   // this in principle should stay on a T0 PID QA, but are just the data prepared for TOF use
2084   TH1F *hStartTimeAT0 = new TH1F("hStartTimeA_T0","StartTime from T0A [ps]",1000,-1000,1000);
2085   fListQAtof->Add(hStartTimeAT0);
2086   TH1F *hStartTimeCT0 = new TH1F("hStartTimeC_T0","StartTime from T0C [ps]",1000,-1000,1000);
2087   fListQAtof->Add(hStartTimeCT0);
2088   TH1F *hStartTimeACT0 = new TH1F("hStartTimeAC_T0","StartTime from T0AC [ps]",1000,-1000,1000);;
2089   fListQAtof->Add(hStartTimeACT0);
2090 }
2091
2092
2093 //______________________________________________________________________________
2094 void AliAnalysisTaskPIDqa::SetupT0qa()
2095 {
2096   //
2097   // Create the T0 qa objects
2098   //
2099
2100   // these are similar to plots inside TOFqa, but these are for all events
2101   TH1F *hStartTimeAT0 = new TH1F("hStartTimeA_T0","StartTime from T0A [ps]",1000,-1000,1000);
2102   fListQAt0->Add(hStartTimeAT0);
2103   TH1F *hStartTimeCT0 = new TH1F("hStartTimeC_T0","StartTime from T0C [ps]",1000,-1000,1000);
2104   fListQAt0->Add(hStartTimeCT0);
2105   TH1F *hStartTimeACT0 = new TH1F("hStartTimeAC_T0","StartTime from T0AC [ps]",1000,-1000,1000);;
2106   fListQAt0->Add(hStartTimeACT0);
2107
2108   TH1F *hnTracksAtT0 = new TH1F("hnTracksAt_T0","Tracks for events selected for T0",100,0,100);
2109   fListQAt0->Add(hnTracksAtT0);
2110   TH1F *hT0AEff = new TH1F("hT0AEff","Events with T0A vs nTracks",100,0,100);
2111   fListQAt0->Add(hT0AEff);
2112   TH1F *hT0CEff = new TH1F("hT0CEff","Events with T0C vs nTracks",100,0,100);
2113   fListQAt0->Add(hT0CEff);
2114   TH1F *hT0AndEff = new TH1F("hT0AndEff","Events with T0AC (AND) vs nTracks",100,0,100);
2115   fListQAt0->Add(hT0AndEff);
2116   TH1F *hT0OrEff = new TH1F("hT0OrEff","Events with T0AC (OR) vs nTracks",100,0,100);
2117   fListQAt0->Add(hT0OrEff);
2118
2119
2120 }
2121
2122 //______________________________________________________________________________
2123 void AliAnalysisTaskPIDqa::SetupEMCALqa()
2124 {
2125   //
2126   // Create the EMCAL qa objects
2127   //
2128
2129   TVectorD *vX=MakeLogBinning(200,.1,30);
2130
2131   TH2F *hNsigmaPt = new TH2F(Form("hNsigmaPt_EMCAL_%s",AliPID::ParticleName(0)),
2132                              Form("EMCAL n#sigma %s vs. p_{T};p_{T} [GeV]; n#sigma",AliPID::ParticleName(0)),
2133                              vX->GetNrows()-1,vX->GetMatrixArray(),
2134                              200,-10,10);
2135   fListQAemcal->Add(hNsigmaPt);
2136
2137
2138   TH2F *hSigPtEle = new TH2F("hSigPt_EMCAL_Ele",
2139                         "EMCAL signal (E/p) vs. p_{T} for electrons;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
2140                         vX->GetNrows()-1,vX->GetMatrixArray(),
2141                         200,0,2);
2142   fListQAemcal->Add(hSigPtEle);
2143
2144   TH2F *hSigPtPions = new TH2F("hSigPt_EMCAL_Pions",
2145                         "EMCAL signal (E/p) vs. p_{T} for pions;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
2146                         vX->GetNrows()-1,vX->GetMatrixArray(),
2147                         200,0,2);
2148   fListQAemcal->Add(hSigPtPions);
2149
2150   TH2F *hSigPtProtons = new TH2F("hSigPt_EMCAL_Protons",
2151                         "EMCAL signal (E/p) vs. p_{T} for protons;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
2152                         vX->GetNrows()-1,vX->GetMatrixArray(),
2153                         200,0,2);
2154   fListQAemcal->Add(hSigPtProtons);
2155
2156   TH2F *hSigPtAntiProtons = new TH2F("hSigPt_EMCAL_Antiprotons",
2157                         "EMCAL signal (E/p) vs. p_{T} for antiprotons;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
2158                         vX->GetNrows()-1,vX->GetMatrixArray(),
2159                         200,0,2);
2160   fListQAemcal->Add(hSigPtAntiProtons);
2161
2162   delete vX;
2163 }
2164
2165 //______________________________________________________________________________
2166 void AliAnalysisTaskPIDqa::SetupHMPIDqa()
2167 {
2168   //
2169   // Create the HMPID qa objects
2170   //
2171
2172   TVectorD *vX=MakeLogBinning(200,.1,30);
2173
2174   // nSigmas
2175   Int_t nhists=0;
2176   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
2177     if (ispecie==AliPID::kElectron || ispecie==AliPID::kMuon) continue;
2178     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_HMPID_%s",AliPID::ParticleName(ispecie)),
2179                               Form("HMPID n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
2180                               vX->GetNrows()-1,vX->GetMatrixArray(),
2181                               200,-10,10);
2182     fListQAhmpid->AddAt(hNsigmaP, nhists);
2183     ++nhists;
2184   }
2185
2186   // cherenkov angle
2187   TH2F *hCkovAnglevsMom   = new TH2F("hCkovAnglevsMom",  "Cherenkov angle vs momentum",
2188                                      vX->GetNrows()-1,vX->GetMatrixArray(),
2189                                      500,0,1);
2190   fListQAhmpid->AddAt(hCkovAnglevsMom,nhists);
2191
2192   delete vX;
2193 }
2194
2195 //______________________________________________________________________________
2196 void AliAnalysisTaskPIDqa::SetupTOFHMPIDqa()
2197 {
2198   //
2199   // Create the HMPID qa objects
2200   //
2201
2202   TH2F *hCkovAnglevsMomPion   = new TH2F("hCkovAnglevsMom_pion",  "Cherenkov angle vs momentum for pions",500,0,5.,500,0,1);
2203   fListQAtofhmpid->Add(hCkovAnglevsMomPion);
2204
2205   TH2F *hCkovAnglevsMomKaon   = new TH2F("hCkovAnglevsMom_kaon",  "Cherenkov angle vs momentum for kaons",500,0,5.,500,0,1);
2206   fListQAtofhmpid->Add(hCkovAnglevsMomKaon);
2207
2208   TH2F *hCkovAnglevsMomProton = new TH2F("hCkovAnglevsMom_proton","Cherenkov angle vs momentum for protons",500,0,5.,500,0,1);
2209   fListQAtofhmpid->Add(hCkovAnglevsMomProton);
2210
2211
2212 }
2213
2214 //______________________________________________________________________________
2215 void AliAnalysisTaskPIDqa::SetupTPCTOFqa()
2216 {
2217   //
2218   // Create the qa objects for TPC + TOF combination
2219   //
2220
2221   TVectorD *vX=MakeLogBinning(200,.1,30);
2222
2223   //TPC signals after TOF cut
2224   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
2225     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_TOF_%s",AliPID::ParticleName(ispecie)),
2226                               Form("TPC n#sigma %s vs. p (after TOF 3#sigma cut);p_{TPC} [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
2227                               vX->GetNrows()-1,vX->GetMatrixArray(),
2228                               200,-10,10);
2229     fListQAtpctof->Add(hNsigmaP);
2230   }
2231
2232   //TOF signals after TPC cut
2233   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
2234     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_TPC_%s",AliPID::ParticleName(ispecie)),
2235                               Form("TOF n#sigma %s vs. p (after TPC n#sigma cut);p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
2236                               vX->GetNrows()-1,vX->GetMatrixArray(),
2237                               200,-10,10);
2238     fListQAtpctof->Add(hNsigmaP);
2239   }
2240
2241   //EMCAL signal after TOF and TPC cut
2242   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
2243     TH2F *heopPt = new TH2F(Form("heopPt_TOF_TPC_%s",AliPID::ParticleName(ispecie)),
2244                             Form("EMCAL signal (E/p) %s vs. p_{T};p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",AliPID::ParticleName(ispecie)),
2245                             vX->GetNrows()-1,vX->GetMatrixArray(),
2246                             200,0,2);
2247     fListQAtpctof->Add(heopPt);
2248   }
2249
2250   delete vX;
2251 }
2252 //______________________________________________________________________________
2253 void AliAnalysisTaskPIDqa::SetupV0qa()
2254 {
2255   //
2256   // Create the qa objects for V0 Kine cuts
2257   //
2258
2259   TH2F *hArmenteros  = new TH2F("hArmenteros",  "Armenteros plot",200,-1.,1.,200,0.,0.4);
2260   fListQAV0->Add(hArmenteros);
2261
2262 }
2263
2264 //_____________________________________________________________________________
2265 void AliAnalysisTaskPIDqa::SetupQAinfo(){
2266   //
2267   // Setup the info of QA objects
2268   //
2269
2270   TObjArray *arr=new TObjArray;
2271   arr->SetName("TPC_info");
2272   fListQAinfo->Add(arr);
2273 }
2274
2275 //______________________________________________________________________________
2276 TVectorD* AliAnalysisTaskPIDqa::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
2277 {
2278   //
2279   // Make logarithmic binning
2280   // the user has to delete the array afterwards!!!
2281   //
2282
2283   //check limits
2284   if (xmin<1e-20 || xmax<1e-20){
2285     AliError("For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
2286     return MakeLinBinning(nbinsX, xmin, xmax);
2287   }
2288   if (xmax<xmin){
2289     Double_t tmp=xmin;
2290     xmin=xmax;
2291     xmax=tmp;
2292   }
2293   TVectorD *binLim=new TVectorD(nbinsX+1);
2294   Double_t first=xmin;
2295   Double_t last=xmax;
2296   Double_t expMax=TMath::Log(last/first);
2297   for (Int_t i=0; i<nbinsX+1; ++i){
2298     (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
2299   }
2300   return binLim;
2301 }
2302
2303 //______________________________________________________________________________
2304 TVectorD* AliAnalysisTaskPIDqa::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
2305 {
2306   //
2307   // Make linear binning
2308   // the user has to delete the array afterwards!!!
2309   //
2310   if (xmax<xmin){
2311     Double_t tmp=xmin;
2312     xmin=xmax;
2313     xmax=tmp;
2314   }
2315   TVectorD *binLim=new TVectorD(nbinsX+1);
2316   Double_t first=xmin;
2317   Double_t last=xmax;
2318   Double_t binWidth=(last-first)/nbinsX;
2319   for (Int_t i=0; i<nbinsX+1; ++i){
2320     (*binLim)[i]=first+binWidth*(Double_t)i;
2321   }
2322   return binLim;
2323 }
2324
2325 //_____________________________________________________________________________
2326 TVectorD* AliAnalysisTaskPIDqa::MakeArbitraryBinning(const char* bins)
2327 {
2328   //
2329   // Make arbitrary binning, bins separated by a ','
2330   //
2331   TString limits(bins);
2332   if (limits.IsNull()){
2333     AliError("Bin Limit string is empty, cannot add the variable");
2334     return 0x0;
2335   }
2336
2337   TObjArray *arr=limits.Tokenize(",");
2338   Int_t nLimits=arr->GetEntries();
2339   if (nLimits<2){
2340     AliError("Need at leas 2 bin limits, cannot add the variable");
2341     delete arr;
2342     return 0x0;
2343   }
2344
2345   TVectorD *binLimits=new TVectorD(nLimits);
2346   for (Int_t iLim=0; iLim<nLimits; ++iLim){
2347     (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
2348   }
2349
2350   delete arr;
2351   return binLimits;
2352 }
2353