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
  41 #include <AliESDEvent.h>
  42 #include <AliAODEvent.h>
  43 #include <AliESDv0.h>
  44 #include <AliAODv0.h>
  45 #include <AliESDv0KineCuts.h>
  46
  47 #include "AliAnalysisTaskPIDqa.h"
  48
  49
  50 ClassImp(AliAnalysisTaskPIDqa)
  51
  52 //______________________________________________________________________________
  53 AliAnalysisTaskPIDqa::AliAnalysisTaskPIDqa():
  54 AliAnalysisTaskSE(),
  55 fPIDResponse(0x0),
  56 fV0cuts(0x0),
  57 fV0electrons(0x0),
  58 fV0pions(0x0),
  59 fV0kaons(0x0),
  60 fV0protons(0x0),
  61 fListQA(0x0),
  62 fListQAits(0x0),
  63 fListQAitsSA(0x0),
  64 fListQAitsPureSA(0x0),
  65 fListQAtpc(0x0),
  66 fListQAtrd(0x0),
  67 fListQAtof(0x0),
  68 fListQAt0(0x0),
  69 fListQAemcal(0x0),
  70 fListQAhmpid(0x0),
  71 fListQAtofhmpid(0x0),
  72 fListQAtpctof(0x0),
  73 fListQAV0(0x0),
  74 fListQAinfo(0x0)
  75 {
  76   //
  77   // Dummy constructor
  78   //
  79 }
  80
  81 //______________________________________________________________________________
  82 AliAnalysisTaskPIDqa::AliAnalysisTaskPIDqa(const char* name):
  83 AliAnalysisTaskSE(name),
  84 fPIDResponse(0x0),
  85 fV0cuts(0x0),
  86 fV0electrons(0x0),
  87 fV0pions(0x0),
  88 fV0kaons(0x0),
  89 fV0protons(0x0),
  90 fListQA(0x0),
  91 fListQAits(0x0),
  92 fListQAitsSA(0x0),
  93 fListQAitsPureSA(0x0),
  94 fListQAtpc(0x0),
  95 fListQAtrd(0x0),
  96 fListQAtof(0x0),
  97 fListQAt0(0x0),
  98 fListQAemcal(0x0),
  99 fListQAhmpid(0x0),
 100 fListQAtofhmpid(0x0),
 101 fListQAtpctof(0x0),
 102 fListQAV0(0x0),
 103 fListQAinfo(0x0)
 104 {
 105   //
 106   // Default constructor
 107   //
 108   DefineInput(0,TChain::Class());
 109   DefineOutput(1,TList::Class());
 110 }
 111
 112 //______________________________________________________________________________
 113 AliAnalysisTaskPIDqa::~AliAnalysisTaskPIDqa()
 114 {
 115   //
 116   // Destructor
 117   //
 118
 119   delete fV0cuts;
 120   delete fV0electrons;
 121   delete fV0pions;
 122   delete fV0kaons;
 123   delete fV0protons;
 124
 125   if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) delete fListQA;
 126 }
 127
 128 //______________________________________________________________________________
 129 void AliAnalysisTaskPIDqa::UserCreateOutputObjects()
 130 {
 131   //
 132   // Create the output QA objects
 133   //
 134
 135   AliLog::SetClassDebugLevel("AliAnalysisTaskPIDqa",10);
 136
 137   //input hander
 138   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 139   AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
 140   if (!inputHandler) AliFatal("Input handler needed");
 141
 142   //pid response object
 143   fPIDResponse=inputHandler->GetPIDResponse();
 144   if (!fPIDResponse) AliError("PIDResponse object was not created");
 145
 146   // V0 Kine cuts
 147   fV0cuts = new AliESDv0KineCuts;
 148
 149   // V0 PID Obj arrays
 150   fV0electrons = new TObjArray;
 151   fV0pions     = new TObjArray;
 152   fV0kaons     = new TObjArray;
 153   fV0protons   = new TObjArray;
 154
 155   //
 156   fListQA=new TList;
 157   fListQA->SetOwner();
 158
 159   fListQAits=new TList;
 160   fListQAits->SetOwner();
 161   fListQAits->SetName("ITS");
 162
 163   fListQAitsSA=new TList;
 164   fListQAitsSA->SetOwner();
 165   fListQAitsSA->SetName("ITS_SA");
 166
 167   fListQAitsPureSA=new TList;
 168   fListQAitsPureSA->SetOwner();
 169   fListQAitsPureSA->SetName("ITS_PureSA");
 170
 171   fListQAtpc=new TList;
 172   fListQAtpc->SetOwner();
 173   fListQAtpc->SetName("TPC");
 174
 175   fListQAtrd=new TList;
 176   fListQAtrd->SetOwner();
 177   fListQAtrd->SetName("TRD");
 178
 179   fListQAtof=new TList;
 180   fListQAtof->SetOwner();
 181   fListQAtof->SetName("TOF");
 182
 183   fListQAt0=new TList;
 184   fListQAt0->SetOwner();
 185   fListQAt0->SetName("T0");
 186
 187   fListQAemcal=new TList;
 188   fListQAemcal->SetOwner();
 189   fListQAemcal->SetName("EMCAL");
 190
 191   fListQAhmpid=new TList;
 192   fListQAhmpid->SetOwner();
 193   fListQAhmpid->SetName("HMPID");
 194
 195   fListQAtpctof=new TList;
 196   fListQAtpctof->SetOwner();
 197   fListQAtpctof->SetName("TPC_TOF");
 198
 199   fListQAtofhmpid=new TList;
 200   fListQAtofhmpid->SetOwner();
 201   fListQAtofhmpid->SetName("TOF_HMPID");
 202
 203   fListQAV0=new TList;
 204   fListQAV0->SetOwner();
 205   fListQAV0->SetName("V0decay");
 206
 207   fListQAinfo=new TList;
 208   fListQAinfo->SetOwner();
 209   fListQAinfo->SetName("QAinfo");
 210
 211   fListQA->Add(fListQAits);
 212   fListQA->Add(fListQAitsSA);
 213   fListQA->Add(fListQAitsPureSA);
 214   fListQA->Add(fListQAtpc);
 215   fListQA->Add(fListQAtrd);
 216   fListQA->Add(fListQAtof);
 217   fListQA->Add(fListQAt0);
 218   fListQA->Add(fListQAemcal);
 219   fListQA->Add(fListQAhmpid);
 220   fListQA->Add(fListQAtpctof);
 221   fListQA->Add(fListQAtofhmpid);
 222   fListQA->Add(fListQAV0);
 223   fListQA->Add(fListQAinfo);
 224
 225   SetupITSqa();
 226   SetupTPCqa();
 227   SetupTRDqa();
 228   SetupTOFqa();
 229   SetupT0qa();
 230   SetupEMCALqa();
 231   SetupHMPIDqa();
 232   SetupTPCTOFqa();
 233   SetupTOFHMPIDqa();
 234   SetupV0qa();
 235   SetupQAinfo();
 236
 237   PostData(1,fListQA);
 238 }
 239
 240
 241 //______________________________________________________________________________
 242 void AliAnalysisTaskPIDqa::UserExec(Option_t */*option*/)
 243 {
 244   //
 245   // Setup the PID response functions and fill the QA histograms
 246   //
 247
 248   AliVEvent *event=InputEvent();
 249   if (!event||!fPIDResponse) return;
 250
 251   // Start with the V0 task (only possible for ESDs?)
 252   FillV0PIDlist();
 253
 254   FillITSqa();
 255   FillTPCqa();
 256   FillTRDqa();
 257   FillTOFqa();
 258   FillEMCALqa();
 259   FillHMPIDqa();
 260   FillT0qa();
 261
 262   //combined detector QA
 263   FillTPCTOFqa();
 264   FillTOFHMPIDqa();
 265
 266   // Clear the V0 PID arrays
 267   ClearV0PIDlist();
 268
 269   //QA info
 270   FillQAinfo();
 271
 272   PostData(1,fListQA);
 273 }
 274
 275 //______________________________________________________________________________
 276 void  AliAnalysisTaskPIDqa::FillV0PIDlist(){
 277
 278   //
 279   // Fill the PID object arrays holding the pointers to identified particle tracks
 280   //
 281
 282   // Dynamic cast to ESD events (DO NOTHING for AOD events)
 283   AliESDEvent *event = dynamic_cast<AliESDEvent *>(InputEvent());
 284   if ( !event )  return;
 285
 286   if(TString(event->GetBeamType())=="Pb-Pb" || TString(event->GetBeamType())=="A-A"){
 287     fV0cuts->SetMode(AliESDv0KineCuts::kPurity,AliESDv0KineCuts::kPbPb);
 288   }
 289   else{
 290     fV0cuts->SetMode(AliESDv0KineCuts::kPurity,AliESDv0KineCuts::kPP);
 291   }
 292
 293   // V0 selection
 294   // set event
 295   fV0cuts->SetEvent(event);
 296
 297   // loop over V0 particles
 298   for(Int_t iv0=0; iv0<event->GetNumberOfV0s();iv0++){
 299
 300     AliESDv0 *v0 = (AliESDv0 *) event->GetV0(iv0);
 301
 302     if(!v0) continue;
 303     if(v0->GetOnFlyStatus()) continue;
 304
 305     // Get the particle selection
 306     Bool_t foundV0 = kFALSE;
 307     Int_t pdgV0, pdgP, pdgN;
 308
 309     foundV0 = fV0cuts->ProcessV0(v0, pdgV0, pdgP, pdgN);
 310     if(!foundV0) continue;
 311
 312     Int_t iTrackP = v0->GetPindex();  // positive track
 313     Int_t iTrackN = v0->GetNindex();  // negative track
 314
 315     // v0 Armenteros plot (QA)
 316     Float_t armVar[2] = {0.0,0.0};
 317     fV0cuts->Armenteros(v0, armVar);
 318
 319     TH2 *h=(TH2*)fListQAV0->At(0);
 320     if (!h) continue;
 321     h->Fill(armVar[0],armVar[1]);
 322
 323     // fill the Object arrays
 324     // positive particles
 325     if( pdgP == -11){
 326       fV0electrons->Add((AliVTrack*)event->GetTrack(iTrackP));
 327     }
 328     else if( pdgP == 211){
 329       fV0pions->Add((AliVTrack*)event->GetTrack(iTrackP));
 330     }
 331     else if( pdgP == 321){
 332       fV0kaons->Add((AliVTrack*)event->GetTrack(iTrackP));
 333     }
 334     else if( pdgP == 2212){
 335       fV0protons->Add((AliVTrack*)event->GetTrack(iTrackP));
 336     }
 337
 338     // negative particles
 339     if( pdgN == 11){
 340       fV0electrons->Add((AliVTrack*)event->GetTrack(iTrackN));
 341     }
 342     else if( pdgN == -211){
 343       fV0pions->Add((AliVTrack*)event->GetTrack(iTrackN));
 344     }
 345     else if( pdgN == -321){
 346       fV0kaons->Add((AliVTrack*)event->GetTrack(iTrackN));
 347     }
 348     else if( pdgN == -2212){
 349       fV0protons->Add((AliVTrack*)event->GetTrack(iTrackN));
 350     }
 351
 352
 353   }
 354 }
 355 //______________________________________________________________________________
 356 void  AliAnalysisTaskPIDqa::ClearV0PIDlist(){
 357
 358   //
 359   // Clear the PID object arrays
 360   //
 361
 362   fV0electrons->Clear();
 363   fV0pions->Clear();
 364   fV0kaons->Clear();
 365   fV0protons->Clear();
 366
 367 }
 368 //______________________________________________________________________________
 369 void AliAnalysisTaskPIDqa::FillITSqa()
 370 {
 371   //
 372   // Fill PID qa histograms for the ITS
 373   //
 374
 375   AliVEvent *event=InputEvent();
 376
 377   Int_t ntracks=event->GetNumberOfTracks();
 378   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 379     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 380     ULong_t status=track->GetStatus();
 381     // not that nice. status bits not in virtual interface
 382     // ITS refit + ITS pid selection
 383     if (!( ( (status & AliVTrack::kITSrefit)==AliVTrack::kITSrefit ) ||
 384            ! ( (status & AliVTrack::kITSpid  )==AliVTrack::kITSpid   ) )) continue;
 385     Double_t mom=track->P();
 386
 387     TList *theList = 0x0;
 388     if(( (status & AliVTrack::kTPCin)==AliVTrack::kTPCin )){
 389       //ITS+TPC tracks
 390       theList=fListQAits;
 391     }else{
 392       if(!( (status & AliVTrack::kITSpureSA)==AliVTrack::kITSpureSA )){
 393         //ITS Standalone tracks
 394         theList=fListQAitsSA;
 395       }else{
 396         //ITS Pure Standalone tracks
 397         theList=fListQAitsPureSA;
 398       }
 399     }
 400
 401
 402     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
 403       TH2 *h=(TH2*)theList->At(ispecie);
 404       if (!h) continue;
 405       Double_t nSigma=fPIDResponse->NumberOfSigmasITS(track, (AliPID::EParticleType)ispecie);
 406       h->Fill(mom,nSigma);
 407     }
 408     TH2 *h=(TH2*)theList->At(AliPID::kSPECIESC);
 409     if (h) {
 410       Double_t sig=track->GetITSsignal();
 411       h->Fill(mom,sig);
 412     }
 413   }
 414 }
 415
 416 //______________________________________________________________________________
 417 void AliAnalysisTaskPIDqa::FillTPCqa()
 418 {
 419   //
 420   // Fill PID qa histograms for the TPC
 421   //
 422
 423   AliVEvent *event=InputEvent();
 424
 425   Int_t ntracks=event->GetNumberOfTracks();
 426   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 427     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 428
 429     //
 430     //basic track cuts
 431     //
 432     ULong_t status=track->GetStatus();
 433     // not that nice. status bits not in virtual interface
 434     // TPC refit + ITS refit + TPC pid
 435     if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 436         !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
 437
 438     // The TPC pid cut removes the light nuclei (>5 sigma from proton line)
 439     //||        !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  )
 440     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 441     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 442     if (track->GetTPCNclsF()>0) {
 443       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 444     }
 445
 446     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 447
 448     Double_t mom=track->GetTPCmomentum();
 449     // the default scenario
 450     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
 451       TH2 *h=(TH2*)fListQAtpc->At(ispecie);
 452       if (!h) continue;
 453       Double_t nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
 454       h->Fill(mom,nSigma);
 455     }
 456     // the "hybrid" scenario
 457     if (track->GetTPCdEdxInfo()){
 458       for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
 459         TH2 *h=(TH2*)fListQAtpc->At(ispecie+AliPID::kSPECIESC);
 460         if (!h) continue;
 461         Double_t nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie, AliTPCPIDResponse::kdEdxHybrid);
 462         h->Fill(mom,nSigma);
 463       }
 464
 465       // the "OROC" scenario
 466       for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
 467         TH2 *h=(TH2*)fListQAtpc->At(ispecie+2*AliPID::kSPECIESC);
 468         if (!h) continue;
 469         Double_t nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie, AliTPCPIDResponse::kdEdxOROC);
 470         //TSpline3* spline = fPIDResponse->GetTPCResponse().GetCurrentResponseFunction();
 471         //std::cout<<ispecie<<" "<<nSigma<<" phi:"<<track->Phi()<<". "<<std::endl;
 472         //if (spline) {cout<<spline->GetName()<<endl;}
 473         //else {cout<<"NULL spline"<<endl;}
 474         h->Fill(mom,nSigma);
 475         }
 476     }
 477
 478     TH2 *h=(TH2*)fListQAtpc->At(3*AliPID::kSPECIESC);
 479
 480     if (h) {
 481       Double_t sig=track->GetTPCsignal();
 482       h->Fill(mom,sig);
 483     }
 484   }
 485 }
 486
 487 //______________________________________________________________________________
 488 void AliAnalysisTaskPIDqa::FillTRDqa()
 489 {
 490   //
 491   // Fill PID qa histograms for the TRD
 492   //
 493   AliVEvent *event=InputEvent();
 494   Int_t ntracks = event->GetNumberOfTracks();
 495   for(Int_t itrack = 0; itrack <  ntracks; itrack++){
 496     AliVTrack *track = (AliVTrack *)event->GetTrack(itrack);
 497
 498     //
 499     //basic track cuts
 500     //
 501     ULong_t status=track->GetStatus();
 502     // not that nice. status bits not in virtual interface
 503     // TPC refit + ITS refit + TPC pid + TRD out
 504     if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 505         !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
 506 //         !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  ) || //removes light nuclei. So it is out for the moment
 507         !( (status & AliVTrack::kTRDout  ) == AliVTrack::kTRDout  )) continue;
 508
 509     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 510     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 511     if (track->GetTPCNclsF()>0) {
 512       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 513     }
 514
 515     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 516
 517     Double_t likelihoods[AliPID::kSPECIES];
 518     if(fPIDResponse->ComputeTRDProbability(track, AliPID::kSPECIES, likelihoods) != AliPIDResponse::kDetPidOk) continue;
 519     Int_t ntracklets = 0;
 520     Double_t momentum = -1.;
 521     for(Int_t itl = 0; itl < 6; itl++)
 522       if(track->GetTRDmomentum(itl) > 0.){
 523         ntracklets++;
 524         if(momentum < 0) momentum = track->GetTRDmomentum(itl);
 525     }
 526     for(Int_t ispecie = 0; ispecie < AliPID::kSPECIES; ispecie++){
 527       TH2F *hLike = (TH2F *)fListQAtrd->At(ntracklets*AliPID::kSPECIES+ispecie);
 528       if (hLike) hLike->Fill(momentum,likelihoods[ispecie]);
 529     }
 530   }
 531 }
 532
 533 //______________________________________________________________________________
 534 void AliAnalysisTaskPIDqa::FillTOFqa()
 535 {
 536   //
 537   // Fill TOF information
 538   //
 539   AliVEvent *event=InputEvent();
 540
 541   Int_t ntracks=event->GetNumberOfTracks();
 542   Int_t tracksAtTof = 0;
 543   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 544     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 545
 546     //
 547     //basic track cuts
 548     //
 549     ULong_t status=track->GetStatus();
 550     // TPC refit + ITS refit +
 551     // TOF out + kTIME
 552     // kTIME
 553     // (we don't use kTOFmismatch because it depends on TPC and kTOFpid because it prevents light nuclei
 554     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 555         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
 556         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||
 557         //        !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||
 558         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;
 559
 560     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 561     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 562     if (track->GetTPCNclsF()>0) {
 563       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 564     }
 565
 566     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 567
 568     tracksAtTof++;
 569
 570     Double_t mom=track->P();
 571
 572     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
 573       TH2 *h=(TH2*)fListQAtof->At(ispecie);
 574       if (!h) continue;
 575       Double_t nSigma=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
 576       h->Fill(mom,nSigma);
 577     }
 578
 579     TH2 *h=(TH2*)fListQAtof->FindObject("hSigP_TOF");
 580     if (h) {
 581       Double_t sig=track->GetTOFsignal()/1000.;
 582       h->Fill(mom,sig);
 583     }
 584
 585     Int_t mask = fPIDResponse->GetTOFResponse().GetStartTimeMask(mom);
 586     ((TH1F*)fListQAtof->FindObject("hStartTimeMask_TOF"))->Fill((Double_t)(mask+0.5));
 587
 588     if (mom >= 0.75 && mom <= 1.25 ) {
 589       Double_t nsigma= fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)AliPID::kPion);
 590       if (mask == 0) {
 591         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-Fill"))->Fill(nsigma);
 592       } else if (mask == 1) {
 593         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-TOF"))->Fill(nsigma);
 594       } else if ( (mask == 2) || (mask == 4) || (mask == 6) ) {
 595         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-T0"))->Fill(nsigma);
 596       } else {
 597         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Pion_T0-Best"))->Fill(nsigma);
 598       }
 599       if (mask & 0x1) { //at least TOF-T0 present
 600         Double_t delta=0;
 601         (void)fPIDResponse->GetSignalDelta((AliPIDResponse::EDetector)AliPIDResponse::kTOF,track,(AliPID::EParticleType)AliPID::kPion,delta);
 602         ((TH1F*)fListQAtof->FindObject("hDelta_TOF_Pion"))->Fill(delta);
 603       }
 604     }
 605
 606     Double_t res = (Double_t)fPIDResponse->GetTOFResponse().GetStartTimeRes(mom);
 607     ((TH1F*)fListQAtof->FindObject("hStartTimeRes_TOF"))->Fill(res);
 608
 609     Double_t startTimeT0 = event->GetT0TOF(0);
 610     if (startTimeT0 < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeAC_T0"))->Fill(startTimeT0);
 611     else {
 612       startTimeT0 = event->GetT0TOF(1);
 613       if (startTimeT0 < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeA_T0"))->Fill(startTimeT0);
 614       startTimeT0 = event->GetT0TOF(2);
 615       if (startTimeT0 < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeC_T0"))->Fill(startTimeT0);
 616     }
 617   }
 618   if (tracksAtTof > 0) {
 619     ((TH1F* )fListQAtof->FindObject("hnTracksAt_TOF"))->Fill(tracksAtTof);
 620     Int_t mask = fPIDResponse->GetTOFResponse().GetStartTimeMask(5.);
 621     if (mask & 0x1) ((TH1F*)fListQAtof->FindObject("hT0MakerEff"))->Fill(tracksAtTof);
 622   }
 623 }
 624
 625 //______________________________________________________________________________
 626 void AliAnalysisTaskPIDqa::FillT0qa()
 627 {
 628   //
 629   // Fill TOF information
 630   //
 631   AliVEvent *event=InputEvent();
 632
 633   Int_t ntracks=event->GetNumberOfTracks();
 634
 635   Int_t tracksAtT0 = 0;
 636
 637   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 638     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 639
 640     //
 641     //basic track cuts
 642     //
 643     ULong_t status=track->GetStatus();
 644     // TPC refit + ITS refit +
 645     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 646         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
 647     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 648     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 649     if (track->GetTPCNclsF()>0) {
 650       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 651     }
 652     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 653
 654     tracksAtT0++;
 655   }
 656
 657   Bool_t t0A = kFALSE;
 658   Bool_t t0C = kFALSE;
 659   Bool_t t0And = kFALSE;
 660   Double_t startTimeT0 = event->GetT0TOF(0);     // AND
 661   if (startTimeT0 < 90000) {
 662     t0And = kTRUE;
 663     ((TH1F*)fListQAt0->FindObject("hStartTimeAC_T0"))->Fill(startTimeT0);
 664     }
 665   startTimeT0 = event->GetT0TOF(1);             // T0A
 666   if (startTimeT0 < 90000) {
 667     t0A = kTRUE;
 668     ((TH1F*)fListQAt0->FindObject("hStartTimeA_T0"))->Fill(startTimeT0);
 669
 670   }
 671   startTimeT0 = event->GetT0TOF(2);             // T0C
 672   if (startTimeT0 < 90000) {
 673     t0C = kTRUE;
 674     ((TH1F*)fListQAt0->FindObject("hStartTimeC_T0"))->Fill(startTimeT0);
 675   }
 676
 677   ((TH1F* )fListQAt0->FindObject("hnTracksAt_T0"))->Fill(tracksAtT0);
 678   if (t0A) ((TH1F*)fListQAt0->FindObject("hT0AEff"))->Fill(tracksAtT0);
 679   if (t0C) ((TH1F*)fListQAt0->FindObject("hT0CEff"))->Fill(tracksAtT0);
 680   if (t0And) ((TH1F*)fListQAt0->FindObject("hT0AndEff"))->Fill(tracksAtT0);
 681   if (t0A || t0C) ((TH1F*)fListQAt0->FindObject("hT0OrEff"))->Fill(tracksAtT0);
 682 }
 683
 684
 685 //______________________________________________________________________________
 686 void AliAnalysisTaskPIDqa::FillEMCALqa()
 687 {
 688   //
 689   // Fill PID qa histograms for the EMCAL
 690   //
 691
 692   AliVEvent *event=InputEvent();
 693
 694   Int_t ntracks=event->GetNumberOfTracks();
 695   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 696     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 697
 698     //
 699     //basic track cuts
 700     //
 701     ULong_t status=track->GetStatus();
 702     // not that nice. status bits not in virtual interface
 703     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
 704
 705     Double_t pt=track->Pt();
 706
 707     //EMCAL nSigma (only for electrons at the moment)
 708     TH2 *h=(TH2*)fListQAemcal->At(0);
 709     if (!h) continue;
 710     Double_t nSigma=fPIDResponse->NumberOfSigmasEMCAL(track, (AliPID::EParticleType)0);
 711     h->Fill(pt,nSigma);
 712
 713   }
 714
 715    //EMCAL signal (E/p vs. pT) for electrons from V0
 716   for(Int_t itrack = 0; itrack < fV0electrons->GetEntries(); itrack++){
 717     AliVTrack *track=(AliVTrack*)fV0electrons->At(itrack);
 718
 719     //
 720     //basic track cuts
 721     //
 722     ULong_t status=track->GetStatus();
 723     // not that nice. status bits not in virtual interface
 724     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
 725
 726     Double_t pt=track->Pt();
 727
 728     TH2 *h=(TH2*)fListQAemcal->At(1);
 729     if (h) {
 730
 731       Int_t nMatchClus = track->GetEMCALcluster();
 732       Double_t mom     = track->P();
 733       Double_t eop     = -1.;
 734
 735       if(nMatchClus > -1){
 736
 737         AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
 738
 739         if(matchedClus){
 740
 741           // matched cluster is EMCAL
 742           if(matchedClus->IsEMCAL()){
 743
 744             Double_t fClsE       = matchedClus->E();
 745             eop                  = fClsE/mom;
 746
 747             h->Fill(pt,eop);
 748
 749           }
 750         }
 751       }
 752     }
 753   }
 754
 755    //EMCAL signal (E/p vs. pT) for pions from V0
 756   for(Int_t itrack = 0; itrack < fV0pions->GetEntries(); itrack++){
 757     AliVTrack *track=(AliVTrack*)fV0pions->At(itrack);
 758
 759     //
 760     //basic track cuts
 761     //
 762     ULong_t status=track->GetStatus();
 763     // not that nice. status bits not in virtual interface
 764     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
 765
 766     Double_t pt=track->Pt();
 767
 768     TH2 *h=(TH2*)fListQAemcal->At(2);
 769     if (h) {
 770
 771       Int_t nMatchClus = track->GetEMCALcluster();
 772       Double_t mom     = track->P();
 773       Double_t eop     = -1.;
 774
 775       if(nMatchClus > -1){
 776
 777         AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
 778
 779         if(matchedClus){
 780
 781           // matched cluster is EMCAL
 782           if(matchedClus->IsEMCAL()){
 783
 784             Double_t fClsE       = matchedClus->E();
 785             eop                  = fClsE/mom;
 786
 787             h->Fill(pt,eop);
 788
 789           }
 790         }
 791       }
 792     }
 793   }
 794
 795    //EMCAL signal (E/p vs. pT) for protons from V0
 796   for(Int_t itrack = 0; itrack < fV0protons->GetEntries(); itrack++){
 797     AliVTrack *track=(AliVTrack*)fV0protons->At(itrack);
 798
 799     //
 800     //basic track cuts
 801     //
 802     ULong_t status=track->GetStatus();
 803     // not that nice. status bits not in virtual interface
 804     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
 805
 806     Double_t pt=track->Pt();
 807
 808     TH2 *hP=(TH2*)fListQAemcal->At(3);
 809     TH2 *hAP=(TH2*)fListQAemcal->At(4);
 810     if (hP && hAP) {
 811
 812       Int_t nMatchClus = track->GetEMCALcluster();
 813       Double_t mom     = track->P();
 814       Int_t charge     = track->Charge();
 815       Double_t eop     = -1.;
 816
 817       if(nMatchClus > -1){
 818
 819         AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
 820
 821         if(matchedClus){
 822
 823           // matched cluster is EMCAL
 824           if(matchedClus->IsEMCAL()){
 825
 826             Double_t fClsE       = matchedClus->E();
 827             eop                  = fClsE/mom;
 828
 829             if(charge > 0)      hP->Fill(pt,eop);
 830             else if(charge < 0) hAP->Fill(pt,eop);
 831
 832           }
 833         }
 834       }
 835     }
 836   }
 837
 838 }
 839
 840
 841 //______________________________________________________________________________
 842 void AliAnalysisTaskPIDqa::FillHMPIDqa()
 843 {
 844   //
 845   // Fill PID qa histograms for the HMPID
 846   //
 847
 848   AliVEvent *event=InputEvent();
 849
 850   Int_t ntracks=event->GetNumberOfTracks();
 851   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 852     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 853
 854     //
 855     //basic track cuts
 856     //
 857     const ULong_t status=track->GetStatus();
 858     // not that nice. status bits not in virtual interface
 859     // TPC refit + ITS refit +
 860     // TOF out + TOFpid +
 861     // kTIME
 862     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 863         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ) continue;
 864
 865     const Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 866     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 867     if (track->GetTPCNclsF()>0) {
 868       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 869     }
 870
 871     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 872
 873     const Double_t mom = track->P();
 874     const Double_t ckovAngle = track->GetHMPIDsignal();
 875
 876     Int_t nhists=0;
 877     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
 878       if (ispecie==AliPID::kElectron || ispecie==AliPID::kMuon) continue;
 879       TH2 *h=(TH2*)fListQAhmpid->At(nhists);
 880       if (!h) {++nhists; continue;}
 881       const Double_t nSigma=fPIDResponse->NumberOfSigmasHMPID(track, (AliPID::EParticleType)ispecie);
 882       h->Fill(mom,nSigma);
 883       ++nhists;
 884     }
 885
 886     TH1F *hThetavsMom = (TH1F*)fListQAhmpid->At(AliPID::kSPECIESC);
 887
 888     if (hThetavsMom) hThetavsMom->Fill(mom,ckovAngle);
 889
 890   }
 891 }
 892 //______________________________________________________________________________
 893 void AliAnalysisTaskPIDqa::FillTOFHMPIDqa()
 894 {
 895   //
 896   // Fill PID qa histograms for the HMPID
 897   //
 898
 899   AliVEvent *event=InputEvent();
 900
 901   Int_t ntracks=event->GetNumberOfTracks();
 902   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 903     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 904
 905     //
 906     //basic track cuts
 907     //
 908     ULong_t status=track->GetStatus();
 909     // not that nice. status bits not in virtual interface
 910     // TPC refit + ITS refit +
 911     // TOF out + TOFpid +
 912     // kTIME
 913     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 914         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
 915         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||
 916         !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||
 917         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;
 918
 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     Double_t mom = track->P();
 928     Double_t ckovAngle = track->GetHMPIDsignal();
 929
 930     Double_t nSigmaTOF[3];
 931     TH1F *h[3];
 932
 933     for (Int_t ispecie=2; ispecie<5; ++ispecie){
 934       //TOF nSigma
 935       nSigmaTOF[ispecie-2]=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
 936       h[ispecie-2] = (TH1F*)fListQAtofhmpid->At(ispecie-2);}
 937
 938     if(TMath::Abs(nSigmaTOF[0])<2)                                                              h[0]->Fill(mom,ckovAngle);
 939
 940     if(TMath::Abs(nSigmaTOF[1])<2 && TMath::Abs(nSigmaTOF[0])>3)                                h[1]->Fill(mom,ckovAngle);
 941
 942     if(TMath::Abs(nSigmaTOF[2])<2 && TMath::Abs(nSigmaTOF[1])>3 && TMath::Abs(nSigmaTOF[0])>3)  h[2]->Fill(mom,ckovAngle);
 943
 944   }
 945
 946 }
 947
 948 //______________________________________________________________________________
 949 void AliAnalysisTaskPIDqa::FillTPCTOFqa()
 950 {
 951   //
 952   // Fill PID qa histograms for the TOF
 953   //   Here also the TPC histograms after TOF selection are filled
 954   //
 955
 956   AliVEvent *event=InputEvent();
 957
 958   Int_t ntracks=event->GetNumberOfTracks();
 959   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 960     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 961
 962     //
 963     //basic track cuts
 964     //
 965     ULong_t status=track->GetStatus();
 966     // not that nice. status bits not in virtual interface
 967     // TPC refit + ITS refit +
 968     // TOF out + TOFpid +
 969     // kTIME
 970     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 971         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
 972 //         !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid ) || //removes light nuclei, so it is out for the moment
 973         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||
 974         !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||
 975         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;
 976
 977     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 978     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 979     if (track->GetTPCNclsF()>0) {
 980       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 981     }
 982
 983     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 984
 985
 986     Double_t mom=track->P();
 987     Double_t momTPC=track->GetTPCmomentum();
 988
 989     for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
 990       //TOF nSigma
 991       Double_t nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
 992       Double_t nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
 993
 994       //TPC after TOF cut
 995       TH2 *h=(TH2*)fListQAtpctof->At(ispecie);
 996       if (h && TMath::Abs(nSigmaTOF)<3.) h->Fill(momTPC,nSigmaTPC);
 997
 998       //TOF after TPC cut
 999       h=(TH2*)fListQAtpctof->At(ispecie+AliPID::kSPECIESC);
1000       if (h && TMath::Abs(nSigmaTPC)<3.) h->Fill(mom,nSigmaTOF);
1001
1002       //EMCAL after TOF and TPC cut
1003       h=(TH2*)fListQAtpctof->At(ispecie+2*AliPID::kSPECIESC);
1004       if (h && TMath::Abs(nSigmaTOF)<3. && TMath::Abs(nSigmaTPC)<3. ){
1005
1006         Int_t nMatchClus = track->GetEMCALcluster();
1007         Double_t pt      = track->Pt();
1008         Double_t eop     = -1.;
1009
1010         if(nMatchClus > -1){
1011
1012           AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
1013
1014           if(matchedClus){
1015
1016             // matched cluster is EMCAL
1017             if(matchedClus->IsEMCAL()){
1018
1019               Double_t fClsE       = matchedClus->E();
1020               eop                  = fClsE/mom;
1021
1022               h->Fill(pt,eop);
1023
1024
1025             }
1026           }
1027         }
1028       }
1029     }
1030   }
1031 }
1032
1033 //_____________________________________________________________________________
1034 void AliAnalysisTaskPIDqa::FillQAinfo()
1035 {
1036   //
1037   // Fill the QA information
1038   //
1039
1040
1041   //TPC QA info
1042   TObjArray *arrTPC=static_cast<TObjArray*>(fListQAinfo->At(0));
1043   if (fPIDResponse && arrTPC){
1044     AliTPCPIDResponse &tpcResp=fPIDResponse->GetTPCResponse();
1045     // fill spline names
1046     if (!arrTPC->UncheckedAt(0)){
1047
1048       TObjArray *arrTPCsplineNames=new TObjArray(AliPID::kSPECIESC);
1049       arrTPCsplineNames->SetOwner();
1050       arrTPCsplineNames->SetName("TPC_spline_names");
1051       arrTPC->AddAt(arrTPCsplineNames,0);
1052
1053       for (Int_t iresp=0; iresp<AliPID::kSPECIESC; ++iresp){
1054         const TObject *o=tpcResp.GetResponseFunction((AliPID::EParticleType)iresp);
1055         if (!o) continue;
1056         arrTPCsplineNames->Add(new TObjString(Form("%02d: %s",iresp, o->GetName())));
1057       }
1058     }
1059
1060     // tpc response config
1061     if (!arrTPC->UncheckedAt(1)){
1062
1063       TObjArray *arrTPCconfigInfo=new TObjArray;
1064       arrTPCconfigInfo->SetOwner();
1065       arrTPCconfigInfo->SetName("TPC_config_info");
1066       arrTPC->AddAt(arrTPCconfigInfo,1);
1067
1068       TObjString *ostr=0x0;
1069       ostr=new TObjString;
1070       ostr->String().Form("Eta Corr map: %s", tpcResp.GetEtaCorrMap()?tpcResp.GetEtaCorrMap()->GetName():"none");
1071       arrTPCconfigInfo->Add(ostr);
1072
1073       ostr=new TObjString;
1074       ostr->String().Form("Sigma Par map: %s", tpcResp.GetSigmaPar1Map()?tpcResp.GetSigmaPar1Map()->GetName():"none");
1075       arrTPCconfigInfo->Add(ostr);
1076
1077       ostr=new TObjString;
1078       ostr->String().Form("MIP: %.2f", tpcResp.GetMIP());
1079       arrTPCconfigInfo->Add(ostr);
1080
1081       ostr=new TObjString;
1082       ostr->String().Form("Res: Def %.3g (%.3g) : AllHigh %.3g (%.3g) : OROC high %.3g (%.3g)",
1083                           tpcResp.GetRes0(AliTPCPIDResponse::kDefault), tpcResp.GetResN2(AliTPCPIDResponse::kDefault),
1084                           tpcResp.GetRes0(AliTPCPIDResponse::kALLhigh), tpcResp.GetResN2(AliTPCPIDResponse::kALLhigh),
1085                           tpcResp.GetRes0(AliTPCPIDResponse::kOROChigh), tpcResp.GetResN2(AliTPCPIDResponse::kOROChigh)
1086                          );
1087       arrTPCconfigInfo->Add(ostr);
1088     }
1089   }
1090 }
1091
1092 //______________________________________________________________________________
1093 void AliAnalysisTaskPIDqa::SetupITSqa()
1094 {
1095   //
1096   // Create the ITS qa objects
1097   //
1098
1099   TVectorD *vX=MakeLogBinning(200,.1,30);
1100
1101   //ITS+TPC tracks
1102   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1103     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_ITS_%s",AliPID::ParticleName(ispecie)),
1104                               Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1105                               vX->GetNrows()-1,vX->GetMatrixArray(),
1106                               200,-10,10);
1107     fListQAits->Add(hNsigmaP);
1108   }
1109   TH2F *hSig = new TH2F("hSigP_ITS",
1110                         "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
1111                         vX->GetNrows()-1,vX->GetMatrixArray(),
1112                         300,0,300);
1113   fListQAits->Add(hSig);
1114
1115   //ITS Standalone tracks
1116   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1117     TH2F *hNsigmaPSA = new TH2F(Form("hNsigmaP_ITSSA_%s",AliPID::ParticleName(ispecie)),
1118                                 Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1119                                 vX->GetNrows()-1,vX->GetMatrixArray(),
1120                                 200,-10,10);
1121     fListQAitsSA->Add(hNsigmaPSA);
1122   }
1123   TH2F *hSigSA = new TH2F("hSigP_ITSSA",
1124                           "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
1125                           vX->GetNrows()-1,vX->GetMatrixArray(),
1126                           300,0,300);
1127   fListQAitsSA->Add(hSigSA);
1128
1129   //ITS Pure Standalone tracks
1130   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1131     TH2F *hNsigmaPPureSA = new TH2F(Form("hNsigmaP_ITSPureSA_%s",AliPID::ParticleName(ispecie)),
1132                                     Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1133                                     vX->GetNrows()-1,vX->GetMatrixArray(),
1134                                     200,-10,10);
1135     fListQAitsPureSA->Add(hNsigmaPPureSA);
1136   }
1137   TH2F *hSigPureSA = new TH2F("hSigP_ITSPureSA",
1138                               "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
1139                               vX->GetNrows()-1,vX->GetMatrixArray(),
1140                               300,0,300);
1141   fListQAitsPureSA->Add(hSigPureSA);
1142
1143   delete vX;
1144 }
1145
1146 //______________________________________________________________________________
1147 void AliAnalysisTaskPIDqa::SetupTPCqa()
1148 {
1149   //
1150   // Create the TPC qa objects
1151   //
1152
1153   TVectorD *vX=MakeLogBinning(200,.1,30);
1154
1155   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1156     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_%s",AliPID::ParticleName(ispecie)),
1157                               Form("TPC n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1158                               vX->GetNrows()-1,vX->GetMatrixArray(),
1159                               200,-10,10);
1160     fListQAtpc->Add(hNsigmaP);
1161   }
1162
1163   // the "hybrid" scenario
1164   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1165     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_%s_Hybrid",AliPID::ParticleName(ispecie)),
1166                               Form("TPC n#sigma %s vs. p (Hybrid gain scenario);p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1167                               vX->GetNrows()-1,vX->GetMatrixArray(),
1168                               200,-10,10);
1169     fListQAtpc->Add(hNsigmaP);
1170   }
1171
1172   // the "OROC high" scenario
1173   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1174     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_%s_OROChigh",AliPID::ParticleName(ispecie)),
1175                               Form("TPC n#sigma %s vs. p (OROChigh gain scenario);p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1176                               vX->GetNrows()-1,vX->GetMatrixArray(),
1177                               200,-10,10);
1178     fListQAtpc->Add(hNsigmaP);
1179   }
1180
1181
1182
1183   TH2F *hSig = new TH2F("hSigP_TPC",
1184                         "TPC signal vs. p;p [GeV]; TPC signal [arb. units]",
1185                         vX->GetNrows()-1,vX->GetMatrixArray(),
1186                         500,0,1000);
1187   fListQAtpc->Add(hSig); //3*AliPID::kSPECIESC
1188
1189   delete vX;
1190 }
1191
1192 //______________________________________________________________________________
1193 void AliAnalysisTaskPIDqa::SetupTRDqa()
1194 {
1195   //
1196   // Create the TRD qa objects
1197   //
1198   TVectorD *vX=MakeLogBinning(200,.1,30);
1199   for(Int_t itl = 0; itl < 6; ++itl){
1200     for(Int_t ispecie = 0; ispecie < AliPID::kSPECIES; ispecie++){
1201       TH2F *hLikeP = new TH2F(Form("hLikeP_TRD_%dtls_%s", itl, AliPID::ParticleName(ispecie)),
1202                               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)),
1203                               vX->GetNrows()-1, vX->GetMatrixArray(),
1204                               100, 0., 1.);
1205       fListQAtrd->Add(hLikeP);
1206     }
1207   }
1208   delete vX;
1209 }
1210
1211 //______________________________________________________________________________
1212 void AliAnalysisTaskPIDqa::SetupTOFqa()
1213 {
1214   //
1215   // Create the TOF qa objects
1216   //
1217
1218   TVectorD *vX=MakeLogBinning(200,.1,30);
1219
1220   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1221     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_%s",AliPID::ParticleName(ispecie)),
1222                               Form("TOF n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1223                               vX->GetNrows()-1,vX->GetMatrixArray(),
1224                               200,-10,10);
1225     fListQAtof->Add(hNsigmaP);
1226   }
1227
1228   TH1F *hnSigT0Fill = new TH1F("hNsigma_TOF_Pion_T0-Fill","TOF n#sigma (Pion) T0-FILL [0.75-1.25. GeV/c]",200,-10,10);
1229   fListQAtof->Add(hnSigT0Fill);
1230   TH1F *hnSigT0T0 = new TH1F("hNsigma_TOF_Pion_T0-T0","TOF n#sigma (Pion) T0-T0 [0.75-1.25 GeV/c]",200,-10,10);
1231   fListQAtof->Add(hnSigT0T0);
1232   TH1F *hnSigT0TOF = new TH1F("hNsigma_TOF_Pion_T0-TOF","TOF n#sigma (Pion) T0-TOF [0.75-1.25 GeV/c]",200,-10,10);
1233   fListQAtof->Add(hnSigT0TOF);
1234   TH1F *hnSigT0Best = new TH1F("hNsigma_TOF_Pion_T0-Best","TOF n#sigma (Pion) T0-Best [0.75-1.25 GeV/c]",200,-10,10);
1235   fListQAtof->Add(hnSigT0Best);
1236   TH1F *hnDeltaPi = new TH1F("hDelta_TOF_Pion","DeltaT (Pion) [0.75-1.25 GeV/c]",50,-500,500);
1237   fListQAtof->Add(hnDeltaPi);
1238
1239   TH2F *hSig = new TH2F("hSigP_TOF",
1240                         "TOF signal vs. p;p [GeV]; TOF signal [ns]",
1241                         vX->GetNrows()-1,vX->GetMatrixArray(),
1242                         300,0,30);
1243
1244   delete vX;
1245
1246   fListQAtof->Add(hSig);
1247
1248   TH1F *hStartTimeMaskTOF = new TH1F("hStartTimeMask_TOF","StartTime mask",8,0,8);
1249   fListQAtof->Add(hStartTimeMaskTOF);
1250   TH1F *hStartTimeResTOF = new TH1F("hStartTimeRes_TOF","StartTime resolution [ps]",100,0,500);
1251   fListQAtof->Add(hStartTimeResTOF);
1252
1253   TH1F *hnTracksAtTOF = new TH1F("hnTracksAt_TOF","Matched tracks at TOF",100,0,100);
1254   fListQAtof->Add(hnTracksAtTOF);
1255   TH1F *hT0MakerEff = new TH1F("hT0MakerEff","Events with T0-TOF vs nTracks",100,0,100);
1256   fListQAtof->Add(hT0MakerEff);
1257
1258   // this in principle should stay on a T0 PID QA, but are just the data prepared for TOF use
1259   TH1F *hStartTimeAT0 = new TH1F("hStartTimeA_T0","StartTime from T0A [ps]",1000,-1000,1000);
1260   fListQAtof->Add(hStartTimeAT0);
1261   TH1F *hStartTimeCT0 = new TH1F("hStartTimeC_T0","StartTime from T0C [ps]",1000,-1000,1000);
1262   fListQAtof->Add(hStartTimeCT0);
1263   TH1F *hStartTimeACT0 = new TH1F("hStartTimeAC_T0","StartTime from T0AC [ps]",1000,-1000,1000);;
1264   fListQAtof->Add(hStartTimeACT0);
1265 }
1266
1267
1268 //______________________________________________________________________________
1269 void AliAnalysisTaskPIDqa::SetupT0qa()
1270 {
1271   //
1272   // Create the T0 qa objects
1273   //
1274
1275   // these are similar to plots inside TOFqa, but these are for all events
1276   TH1F *hStartTimeAT0 = new TH1F("hStartTimeA_T0","StartTime from T0A [ps]",1000,-1000,1000);
1277   fListQAt0->Add(hStartTimeAT0);
1278   TH1F *hStartTimeCT0 = new TH1F("hStartTimeC_T0","StartTime from T0C [ps]",1000,-1000,1000);
1279   fListQAt0->Add(hStartTimeCT0);
1280   TH1F *hStartTimeACT0 = new TH1F("hStartTimeAC_T0","StartTime from T0AC [ps]",1000,-1000,1000);;
1281   fListQAt0->Add(hStartTimeACT0);
1282
1283   TH1F *hnTracksAtT0 = new TH1F("hnTracksAt_T0","Tracks for events selected for T0",100,0,100);
1284   fListQAt0->Add(hnTracksAtT0);
1285   TH1F *hT0AEff = new TH1F("hT0AEff","Events with T0A vs nTracks",100,0,100);
1286   fListQAt0->Add(hT0AEff);
1287   TH1F *hT0CEff = new TH1F("hT0CEff","Events with T0C vs nTracks",100,0,100);
1288   fListQAt0->Add(hT0CEff);
1289   TH1F *hT0AndEff = new TH1F("hT0AndEff","Events with T0AC (AND) vs nTracks",100,0,100);
1290   fListQAt0->Add(hT0AndEff);
1291   TH1F *hT0OrEff = new TH1F("hT0OrEff","Events with T0AC (OR) vs nTracks",100,0,100);
1292   fListQAt0->Add(hT0OrEff);
1293
1294
1295 }
1296
1297 //______________________________________________________________________________
1298 void AliAnalysisTaskPIDqa::SetupEMCALqa()
1299 {
1300   //
1301   // Create the EMCAL qa objects
1302   //
1303
1304   TVectorD *vX=MakeLogBinning(200,.1,30);
1305
1306   TH2F *hNsigmaPt = new TH2F(Form("hNsigmaPt_EMCAL_%s",AliPID::ParticleName(0)),
1307                              Form("EMCAL n#sigma %s vs. p_{T};p_{T} [GeV]; n#sigma",AliPID::ParticleName(0)),
1308                              vX->GetNrows()-1,vX->GetMatrixArray(),
1309                              200,-10,10);
1310   fListQAemcal->Add(hNsigmaPt);
1311
1312
1313   TH2F *hSigPtEle = new TH2F("hSigPt_EMCAL_Ele",
1314                         "EMCAL signal (E/p) vs. p_{T} for electrons;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
1315                         vX->GetNrows()-1,vX->GetMatrixArray(),
1316                         200,0,2);
1317   fListQAemcal->Add(hSigPtEle);
1318
1319   TH2F *hSigPtPions = new TH2F("hSigPt_EMCAL_Pions",
1320                         "EMCAL signal (E/p) vs. p_{T} for pions;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
1321                         vX->GetNrows()-1,vX->GetMatrixArray(),
1322                         200,0,2);
1323   fListQAemcal->Add(hSigPtPions);
1324
1325   TH2F *hSigPtProtons = new TH2F("hSigPt_EMCAL_Protons",
1326                         "EMCAL signal (E/p) vs. p_{T} for protons;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
1327                         vX->GetNrows()-1,vX->GetMatrixArray(),
1328                         200,0,2);
1329   fListQAemcal->Add(hSigPtProtons);
1330
1331   TH2F *hSigPtAntiProtons = new TH2F("hSigPt_EMCAL_Antiprotons",
1332                         "EMCAL signal (E/p) vs. p_{T} for antiprotons;p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
1333                         vX->GetNrows()-1,vX->GetMatrixArray(),
1334                         200,0,2);
1335   fListQAemcal->Add(hSigPtAntiProtons);
1336
1337   delete vX;
1338 }
1339
1340 //______________________________________________________________________________
1341 void AliAnalysisTaskPIDqa::SetupHMPIDqa()
1342 {
1343   //
1344   // Create the HMPID qa objects
1345   //
1346
1347   TVectorD *vX=MakeLogBinning(200,.1,30);
1348
1349   // nSigmas
1350   Int_t nhists=0;
1351   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
1352     if (ispecie==AliPID::kElectron || ispecie==AliPID::kMuon) continue;
1353     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_HMPID_%s",AliPID::ParticleName(ispecie)),
1354                               Form("HMPID n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1355                               vX->GetNrows()-1,vX->GetMatrixArray(),
1356                               200,-10,10);
1357     fListQAhmpid->AddAt(hNsigmaP, nhists);
1358     ++nhists;
1359   }
1360
1361   // cherenkov angle
1362   TH2F *hCkovAnglevsMom   = new TH2F("hCkovAnglevsMom",  "Cherenkov angle vs momentum",
1363                                      vX->GetNrows()-1,vX->GetMatrixArray(),
1364                                      500,0,1);
1365   fListQAhmpid->AddAt(hCkovAnglevsMom,nhists);
1366
1367   delete vX;
1368 }
1369
1370 //______________________________________________________________________________
1371 void AliAnalysisTaskPIDqa::SetupTOFHMPIDqa()
1372 {
1373   //
1374   // Create the HMPID qa objects
1375   //
1376
1377   TH2F *hCkovAnglevsMomPion   = new TH2F("hCkovAnglevsMom_pion",  "Cherenkov angle vs momentum for pions",500,0,5.,500,0,1);
1378   fListQAtofhmpid->Add(hCkovAnglevsMomPion);
1379
1380   TH2F *hCkovAnglevsMomKaon   = new TH2F("hCkovAnglevsMom_kaon",  "Cherenkov angle vs momentum for kaons",500,0,5.,500,0,1);
1381   fListQAtofhmpid->Add(hCkovAnglevsMomKaon);
1382
1383   TH2F *hCkovAnglevsMomProton = new TH2F("hCkovAnglevsMom_proton","Cherenkov angle vs momentum for protons",500,0,5.,500,0,1);
1384   fListQAtofhmpid->Add(hCkovAnglevsMomProton);
1385
1386
1387 }
1388
1389 //______________________________________________________________________________
1390 void AliAnalysisTaskPIDqa::SetupTPCTOFqa()
1391 {
1392   //
1393   // Create the qa objects for TPC + TOF combination
1394   //
1395
1396   TVectorD *vX=MakeLogBinning(200,.1,30);
1397
1398   //TPC signals after TOF cut
1399   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1400     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_TOF_%s",AliPID::ParticleName(ispecie)),
1401                               Form("TPC n#sigma %s vs. p (after TOF 3#sigma cut);p_{TPC} [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1402                               vX->GetNrows()-1,vX->GetMatrixArray(),
1403                               200,-10,10);
1404     fListQAtpctof->Add(hNsigmaP);
1405   }
1406
1407   //TOF signals after TPC cut
1408   for (Int_t ispecie=0; ispecie<AliPID::kSPECIESC; ++ispecie){
1409     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_TPC_%s",AliPID::ParticleName(ispecie)),
1410                               Form("TOF n#sigma %s vs. p (after TPC n#sigma cut);p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
1411                               vX->GetNrows()-1,vX->GetMatrixArray(),
1412                               200,-10,10);
1413     fListQAtpctof->Add(hNsigmaP);
1414   }
1415
1416   //EMCAL signal after TOF and TPC cut
1417   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
1418     TH2F *heopPt = new TH2F(Form("heopPt_TOF_TPC_%s",AliPID::ParticleName(ispecie)),
1419                             Form("EMCAL signal (E/p) %s vs. p_{T};p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",AliPID::ParticleName(ispecie)),
1420                             vX->GetNrows()-1,vX->GetMatrixArray(),
1421                             200,0,2);
1422     fListQAtpctof->Add(heopPt);
1423   }
1424
1425   delete vX;
1426 }
1427 //______________________________________________________________________________
1428 void AliAnalysisTaskPIDqa::SetupV0qa()
1429 {
1430   //
1431   // Create the qa objects for V0 Kine cuts
1432   //
1433
1434   TH2F *hArmenteros  = new TH2F("hArmenteros",  "Armenteros plot",200,-1.,1.,200,0.,0.4);
1435   fListQAV0->Add(hArmenteros);
1436
1437 }
1438
1439 //_____________________________________________________________________________
1440 void AliAnalysisTaskPIDqa::SetupQAinfo(){
1441   //
1442   // Setup the info of QA objects
1443   //
1444
1445   TObjArray *arr=new TObjArray;
1446   arr->SetName("TPC_info");
1447   fListQAinfo->Add(arr);
1448 }
1449
1450 //______________________________________________________________________________
1451 TVectorD* AliAnalysisTaskPIDqa::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
1452 {
1453   //
1454   // Make logarithmic binning
1455   // the user has to delete the array afterwards!!!
1456   //
1457
1458   //check limits
1459   if (xmin<1e-20 || xmax<1e-20){
1460     AliError("For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
1461     return MakeLinBinning(nbinsX, xmin, xmax);
1462   }
1463   if (xmax<xmin){
1464     Double_t tmp=xmin;
1465     xmin=xmax;
1466     xmax=tmp;
1467   }
1468   TVectorD *binLim=new TVectorD(nbinsX+1);
1469   Double_t first=xmin;
1470   Double_t last=xmax;
1471   Double_t expMax=TMath::Log(last/first);
1472   for (Int_t i=0; i<nbinsX+1; ++i){
1473     (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
1474   }
1475   return binLim;
1476 }
1477
1478 //______________________________________________________________________________
1479 TVectorD* AliAnalysisTaskPIDqa::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
1480 {
1481   //
1482   // Make linear binning
1483   // the user has to delete the array afterwards!!!
1484   //
1485   if (xmax<xmin){
1486     Double_t tmp=xmin;
1487     xmin=xmax;
1488     xmax=tmp;
1489   }
1490   TVectorD *binLim=new TVectorD(nbinsX+1);
1491   Double_t first=xmin;
1492   Double_t last=xmax;
1493   Double_t binWidth=(last-first)/nbinsX;
1494   for (Int_t i=0; i<nbinsX+1; ++i){
1495     (*binLim)[i]=first+binWidth*(Double_t)i;
1496   }
1497   return binLim;
1498 }
1499
1500 //_____________________________________________________________________________
1501 TVectorD* AliAnalysisTaskPIDqa::MakeArbitraryBinning(const char* bins)
1502 {
1503   //
1504   // Make arbitrary binning, bins separated by a ','
1505   //
1506   TString limits(bins);
1507   if (limits.IsNull()){
1508     AliError("Bin Limit string is empty, cannot add the variable");
1509     return 0x0;
1510   }
1511
1512   TObjArray *arr=limits.Tokenize(",");
1513   Int_t nLimits=arr->GetEntries();
1514   if (nLimits<2){
1515     AliError("Need at leas 2 bin limits, cannot add the variable");
1516     delete arr;
1517     return 0x0;
1518   }
1519
1520   TVectorD *binLimits=new TVectorD(nLimits);
1521   for (Int_t iLim=0; iLim<nLimits; ++iLim){
1522     (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
1523   }
1524
1525   delete arr;
1526   return binLimits;
1527 }
1528