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
  43 #include "AliAnalysisTaskPIDqa.h"
  44
  45
  46 ClassImp(AliAnalysisTaskPIDqa)
  47
  48 //______________________________________________________________________________
  49 AliAnalysisTaskPIDqa::AliAnalysisTaskPIDqa():
  50 AliAnalysisTaskSE(),
  51 fPIDResponse(0x0),
  52 fListQA(0x0),
  53 fListQAits(0x0),
  54 fListQAitsSA(0x0),
  55 fListQAitsPureSA(0x0),
  56 fListQAtpc(0x0),
  57 fListQAtrd(0x0),
  58 fListQAtof(0x0),
  59 fListQAemcal(0x0),
  60 fListQAtpctof(0x0)
  61 {
  62   //
  63   // Dummy constructor
  64   //
  65 }
  66
  67 //______________________________________________________________________________
  68 AliAnalysisTaskPIDqa::AliAnalysisTaskPIDqa(const char* name):
  69 AliAnalysisTaskSE(name),
  70 fPIDResponse(0x0),
  71 fListQA(0x0),
  72 fListQAits(0x0),
  73 fListQAitsSA(0x0),
  74 fListQAitsPureSA(0x0),
  75 fListQAtpc(0x0),
  76 fListQAtrd(0x0),
  77 fListQAtof(0x0),
  78 fListQAemcal(0x0),
  79 fListQAtpctof(0x0)
  80 {
  81   //
  82   // Default constructor
  83   //
  84   DefineInput(0,TChain::Class());
  85   DefineOutput(1,TList::Class());
  86 }
  87
  88 //______________________________________________________________________________
  89 AliAnalysisTaskPIDqa::~AliAnalysisTaskPIDqa()
  90 {
  91   //
  92   // Destructor
  93   //
  94
  95 }
  96
  97 //______________________________________________________________________________
  98 void AliAnalysisTaskPIDqa::UserCreateOutputObjects()
  99 {
 100   //
 101   // Create the output QA objects
 102   //
 103
 104   AliLog::SetClassDebugLevel("AliAnalysisTaskPIDqa",10);
 105
 106   //input hander
 107   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
 108   AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
 109   if (!inputHandler) AliFatal("Input handler needed");
 110
 111   //pid response object
 112   fPIDResponse=inputHandler->GetPIDResponse();
 113   if (!fPIDResponse) AliError("PIDResponse object was not created");
 114
 115   //
 116   fListQA=new TList;
 117   fListQA->SetOwner();
 118
 119   fListQAits=new TList;
 120   fListQAits->SetOwner();
 121   fListQAits->SetName("ITS");
 122
 123   fListQAitsSA=new TList;
 124   fListQAitsSA->SetOwner();
 125   fListQAitsSA->SetName("ITS_SA");
 126
 127   fListQAitsPureSA=new TList;
 128   fListQAitsPureSA->SetOwner();
 129   fListQAitsPureSA->SetName("ITS_PureSA");
 130
 131   fListQAtpc=new TList;
 132   fListQAtpc->SetOwner();
 133   fListQAtpc->SetName("TPC");
 134
 135   fListQAtrd=new TList;
 136   fListQAtrd->SetOwner();
 137   fListQAtrd->SetName("TRD");
 138
 139   fListQAtof=new TList;
 140   fListQAtof->SetOwner();
 141   fListQAtof->SetName("TOF");
 142
 143   fListQAemcal=new TList;
 144   fListQAemcal->SetOwner();
 145   fListQAemcal->SetName("EMCAL");
 146
 147   fListQAtpctof=new TList;
 148   fListQAtpctof->SetOwner();
 149   fListQAtpctof->SetName("TPC_TOF");
 150
 151   fListQA->Add(fListQAits);
 152   fListQA->Add(fListQAitsSA);
 153   fListQA->Add(fListQAitsPureSA);
 154   fListQA->Add(fListQAtpc);
 155   fListQA->Add(fListQAtrd);
 156   fListQA->Add(fListQAtof);
 157   fListQA->Add(fListQAemcal);
 158   fListQA->Add(fListQAtpctof);
 159
 160   SetupITSqa();
 161   SetupTPCqa();
 162   SetupTRDqa();
 163   SetupTOFqa();
 164   SetupEMCALqa();
 165   SetupTPCTOFqa();
 166
 167   PostData(1,fListQA);
 168 }
 169
 170
 171 //______________________________________________________________________________
 172 void AliAnalysisTaskPIDqa::UserExec(Option_t */*option*/)
 173 {
 174   //
 175   // Setup the PID response functions and fill the QA histograms
 176   //
 177
 178   AliVEvent *event=InputEvent();
 179   if (!event||!fPIDResponse) return;
 180
 181
 182   FillITSqa();
 183   FillTPCqa();
 184   FillTRDqa();
 185   FillTOFqa();
 186   FillEMCALqa();
 187   FillTPCTOFqa();
 188
 189   PostData(1,fListQA);
 190 }
 191
 192 //______________________________________________________________________________
 193 void AliAnalysisTaskPIDqa::FillITSqa()
 194 {
 195   //
 196   // Fill PID qa histograms for the ITS
 197   //
 198
 199   AliVEvent *event=InputEvent();
 200
 201   Int_t ntracks=event->GetNumberOfTracks();
 202   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 203     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 204     ULong_t status=track->GetStatus();
 205     // not that nice. status bits not in virtual interface
 206     // ITS refit + ITS pid selection
 207     if (!( ( (status & AliVTrack::kITSrefit)==AliVTrack::kITSrefit ) ||
 208            ! ( (status & AliVTrack::kITSpid  )==AliVTrack::kITSpid   ) )) continue;
 209     Double_t mom=track->P();
 210
 211     TList *theList = 0x0;
 212     if(( (status & AliVTrack::kTPCin)==AliVTrack::kTPCin )){
 213       //ITS+TPC tracks
 214       theList=fListQAits;
 215     }else{
 216       if(!( (status & AliVTrack::kITSpureSA)==AliVTrack::kITSpureSA )){
 217         //ITS Standalone tracks
 218         theList=fListQAitsSA;
 219       }else{
 220         //ITS Pure Standalone tracks
 221         theList=fListQAitsPureSA;
 222       }
 223     }
 224
 225
 226     for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 227       TH2 *h=(TH2*)theList->At(ispecie);
 228       if (!h) continue;
 229       Double_t nSigma=fPIDResponse->NumberOfSigmasITS(track, (AliPID::EParticleType)ispecie);
 230       h->Fill(mom,nSigma);
 231     }
 232     TH2 *h=(TH2*)theList->At(AliPID::kSPECIES);
 233     if (h) {
 234       Double_t sig=track->GetITSsignal();
 235       h->Fill(mom,sig);
 236     }
 237   }
 238 }
 239
 240 //______________________________________________________________________________
 241 void AliAnalysisTaskPIDqa::FillTPCqa()
 242 {
 243   //
 244   // Fill PID qa histograms for the TPC
 245   //
 246
 247   AliVEvent *event=InputEvent();
 248
 249   Int_t ntracks=event->GetNumberOfTracks();
 250   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 251     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 252
 253     //
 254     //basic track cuts
 255     //
 256     ULong_t status=track->GetStatus();
 257     // not that nice. status bits not in virtual interface
 258     // TPC refit + ITS refit + TPC pid
 259     if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 260         !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
 261         !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  ) ) continue;
 262
 263     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 264     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 265     if (track->GetTPCNclsF()>0) {
 266       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 267     }
 268
 269     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 270
 271     Double_t mom=track->GetTPCmomentum();
 272
 273     for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 274       TH2 *h=(TH2*)fListQAtpc->At(ispecie);
 275       if (!h) continue;
 276       Double_t nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
 277       h->Fill(mom,nSigma);
 278     }
 279
 280     TH2 *h=(TH2*)fListQAtpc->At(AliPID::kSPECIES);
 281     if (h) {
 282       Double_t sig=track->GetTPCsignal();
 283       h->Fill(mom,sig);
 284     }
 285   }
 286 }
 287
 288 //______________________________________________________________________________
 289 void AliAnalysisTaskPIDqa::FillTRDqa()
 290 {
 291   //
 292   // Fill PID qa histograms for the TRD
 293   //
 294   AliVEvent *event=InputEvent();
 295   Int_t ntracks = event->GetNumberOfTracks();
 296   for(Int_t itrack = 0; itrack <  ntracks; itrack++){
 297     AliVTrack *track = (AliVTrack *)event->GetTrack(itrack);
 298
 299     //
 300     //basic track cuts
 301     //
 302     ULong_t status=track->GetStatus();
 303     // not that nice. status bits not in virtual interface
 304     // TPC refit + ITS refit + TPC pid + TRD out
 305     if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 306         !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
 307         !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  ) ||
 308         !( (status & AliVTrack::kTRDout  ) == AliVTrack::kTRDout  )) continue;
 309
 310     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 311     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 312     if (track->GetTPCNclsF()>0) {
 313       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 314     }
 315
 316     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 317
 318     Double_t likelihoods[AliPID::kSPECIES];
 319     if(fPIDResponse->ComputeTRDProbability(track, AliPID::kSPECIES, likelihoods) != AliPIDResponse::kDetPidOk) continue;
 320     Int_t ntracklets = 0;
 321     Double_t momentum = -1.;
 322     for(Int_t itl = 0; itl < 6; itl++)
 323       if(track->GetTRDmomentum(itl) > 0.){
 324         ntracklets++;
 325         if(momentum < 0) momentum = track->GetTRDmomentum(itl);
 326     }
 327     for(Int_t ispecie = 0; ispecie < AliPID::kSPECIES; ispecie++){
 328       TH2F *hLike = (TH2F *)fListQAtrd->At(ntracklets*AliPID::kSPECIES+ispecie);
 329       if (hLike) hLike->Fill(momentum,likelihoods[ispecie]);
 330     }
 331   }
 332 }
 333
 334 //______________________________________________________________________________
 335 void AliAnalysisTaskPIDqa::FillTOFqa()
 336 {
 337   AliVEvent *event=InputEvent();
 338
 339
 340
 341   Int_t ntracks=event->GetNumberOfTracks();
 342   Int_t tracksAtTof = 0;
 343   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 344     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 345
 346     //
 347     //basic track cuts
 348     //
 349     ULong_t status=track->GetStatus();
 350     // TPC refit + ITS refit +
 351     // TOF out + TOFpid +
 352     // kTIME
 353     // (we don't use kTOFmismatch because it depends on TPC....)
 354     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 355         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
 356         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||
 357         !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||
 358         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;
 359
 360     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 361     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 362     if (track->GetTPCNclsF()>0) {
 363       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 364     }
 365
 366     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 367
 368     tracksAtTof++;
 369
 370     Double_t mom=track->P();
 371
 372     for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 373       TH2 *h=(TH2*)fListQAtof->At(ispecie);
 374       if (!h) continue;
 375       Double_t nSigma=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
 376       h->Fill(mom,nSigma);
 377     }
 378
 379     TH2 *h=(TH2*)fListQAtof->FindObject("hSigP_TOF");
 380     if (h) {
 381       Double_t sig=track->GetTOFsignal();
 382       h->Fill(mom,sig);
 383     }
 384
 385     Double_t mask = (Double_t)fPIDResponse->GetTOFResponse().GetStartTimeMask(mom) + 0.5;
 386     ((TH1F*)fListQAtof->FindObject("hStartTimeMask_TOF"))->Fill(mask);
 387
 388     if (mom >= 1.0 && mom <= 2.0 ) {
 389       Double_t nsigma= fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)AliPID::kKaon);
 390       if (mask == 0) {
 391         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Kaon_T0-Fill"))->Fill(nsigma);
 392       } else if (mask == 1) {
 393         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Kaon_T0-TOF"))->Fill(nsigma);
 394       } else if ( (mask == 2) || (mask == 4) || (mask == 6) ) {
 395         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Kaon_T0-T0"))->Fill(nsigma);
 396       } else {
 397         ((TH1F*)fListQAtof->FindObject("hNsigma_TOF_Kaon_T0-Best"))->Fill(nsigma);
 398       }
 399     }
 400
 401     Double_t res = (Double_t)fPIDResponse->GetTOFResponse().GetStartTimeRes(mom);
 402     ((TH1F*)fListQAtof->FindObject("hStartTimeRes_TOF"))->Fill(res);
 403
 404     AliESDEvent *esd = dynamic_cast<AliESDEvent *>(event);
 405     if (esd) {
 406       Double_t startTime = esd->GetT0TOF(0);
 407       if (startTime < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeAC_T0"))->Fill(startTime);
 408       else {
 409         startTime = esd->GetT0TOF(1);
 410         if (startTime < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeA_T0"))->Fill(startTime);
 411         startTime = esd->GetT0TOF(2);
 412         if (startTime < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeC_T0"))->Fill(startTime);
 413       }
 414     }
 415   }
 416   if (tracksAtTof > 0) {
 417     ((TH1F* )fListQAtof->FindObject("hnTracksAt_TOF"))->Fill(tracksAtTof);
 418     Int_t mask = fPIDResponse->GetTOFResponse().GetStartTimeMask(5.);
 419     if (mask & 0x1) ((TH1F*)fListQAtof->FindObject("hT0MakerEff"))->Fill(tracksAtTof);
 420   }
 421
 422 }
 423
 424
 425 //______________________________________________________________________________
 426 void AliAnalysisTaskPIDqa::FillEMCALqa()
 427 {
 428   //
 429   // Fill PID qa histograms for the EMCAL
 430   //
 431
 432   AliVEvent *event=InputEvent();
 433
 434   Int_t ntracks=event->GetNumberOfTracks();
 435   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 436     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 437
 438     //
 439     //basic track cuts
 440     //
 441     ULong_t status=track->GetStatus();
 442     // not that nice. status bits not in virtual interface
 443     // TPC refit + ITS refit +
 444     // TOF out + TOFpid +
 445     // kTIME
 446     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
 447
 448     Double_t pt=track->Pt();
 449
 450     //EMCAL nSigma (only for electrons at the moment)
 451     TH2 *h=(TH2*)fListQAemcal->At(0);
 452     if (!h) continue;
 453     Double_t nSigma=fPIDResponse->NumberOfSigmasEMCAL(track, (AliPID::EParticleType)0);
 454     h->Fill(pt,nSigma);
 455
 456     //EMCAL signal (E/p vs. pT)
 457     h=(TH2*)fListQAemcal->At(1);
 458     if (h) {
 459
 460       Int_t nMatchClus = track->GetEMCALcluster();
 461       Double_t mom     = track->P();
 462       Double_t eop     = -1.;
 463
 464       if(nMatchClus > -1){
 465
 466         AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
 467
 468         if(matchedClus){
 469
 470           // matched cluster is EMCAL
 471           if(matchedClus->IsEMCAL()){
 472
 473             Double_t fClsE       = matchedClus->E();
 474             eop                  = fClsE/mom;
 475
 476             h->Fill(pt,eop);
 477
 478           }
 479         }
 480       }
 481     }
 482   }
 483 }
 484
 485 //______________________________________________________________________________
 486 void AliAnalysisTaskPIDqa::FillTPCTOFqa()
 487 {
 488   //
 489   // Fill PID qa histograms for the TOF
 490   //   Here also the TPC histograms after TOF selection are filled
 491   //
 492
 493   AliVEvent *event=InputEvent();
 494
 495   Int_t ntracks=event->GetNumberOfTracks();
 496   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 497     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 498
 499     //
 500     //basic track cuts
 501     //
 502     ULong_t status=track->GetStatus();
 503     // not that nice. status bits not in virtual interface
 504     // TPC refit + ITS refit +
 505     // TOF out + TOFpid +
 506     // kTIME
 507     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 508         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
 509         !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid ) ||
 510         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||
 511         !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||
 512         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;
 513
 514     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 515     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 516     if (track->GetTPCNclsF()>0) {
 517       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 518     }
 519
 520     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 521
 522
 523     Double_t mom=track->P();
 524     Double_t momTPC=track->GetTPCmomentum();
 525
 526     for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 527       //TOF nSigma
 528       Double_t nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
 529       Double_t nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
 530
 531       //TPC after TOF cut
 532       TH2 *h=(TH2*)fListQAtpctof->At(ispecie);
 533       if (h && TMath::Abs(nSigmaTOF)<3.) h->Fill(momTPC,nSigmaTPC);
 534
 535       //TOF after TPC cut
 536       h=(TH2*)fListQAtpctof->At(ispecie+AliPID::kSPECIES);
 537       if (h && TMath::Abs(nSigmaTPC)<3.) h->Fill(mom,nSigmaTOF);
 538
 539       //EMCAL after TOF and TPC cut
 540       h=(TH2*)fListQAtpctof->At(ispecie+2*AliPID::kSPECIES);
 541       if (h && TMath::Abs(nSigmaTOF)<3. && TMath::Abs(nSigmaTPC)<3. ){
 542
 543         Int_t nMatchClus = track->GetEMCALcluster();
 544         Double_t pt      = track->Pt();
 545         Double_t eop     = -1.;
 546
 547         if(nMatchClus > -1){
 548
 549           AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
 550
 551           if(matchedClus){
 552
 553             // matched cluster is EMCAL
 554             if(matchedClus->IsEMCAL()){
 555
 556               Double_t fClsE       = matchedClus->E();
 557               eop                  = fClsE/mom;
 558
 559               h->Fill(pt,eop);
 560
 561
 562             }
 563           }
 564         }
 565       }
 566     }
 567   }
 568 }
 569
 570 //______________________________________________________________________________
 571 void AliAnalysisTaskPIDqa::SetupITSqa()
 572 {
 573   //
 574   // Create the ITS qa objects
 575   //
 576
 577   TVectorD *vX=MakeLogBinning(200,.1,30);
 578
 579   //ITS+TPC tracks
 580   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 581     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_ITS_%s",AliPID::ParticleName(ispecie)),
 582                               Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 583                               vX->GetNrows()-1,vX->GetMatrixArray(),
 584                               200,-10,10);
 585     fListQAits->Add(hNsigmaP);
 586   }
 587   TH2F *hSig = new TH2F("hSigP_ITS",
 588                         "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
 589                         vX->GetNrows()-1,vX->GetMatrixArray(),
 590                         300,0,300);
 591   fListQAits->Add(hSig);
 592
 593   //ITS Standalone tracks
 594   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 595     TH2F *hNsigmaPSA = new TH2F(Form("hNsigmaP_ITSSA_%s",AliPID::ParticleName(ispecie)),
 596                                 Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 597                                 vX->GetNrows()-1,vX->GetMatrixArray(),
 598                                 200,-10,10);
 599     fListQAitsSA->Add(hNsigmaPSA);
 600   }
 601   TH2F *hSigSA = new TH2F("hSigP_ITSSA",
 602                           "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
 603                           vX->GetNrows()-1,vX->GetMatrixArray(),
 604                           300,0,300);
 605   fListQAitsSA->Add(hSigSA);
 606
 607   //ITS Pure Standalone tracks
 608   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 609     TH2F *hNsigmaPPureSA = new TH2F(Form("hNsigmaP_ITSPureSA_%s",AliPID::ParticleName(ispecie)),
 610                                     Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 611                                     vX->GetNrows()-1,vX->GetMatrixArray(),
 612                                     200,-10,10);
 613     fListQAitsPureSA->Add(hNsigmaPPureSA);
 614   }
 615   TH2F *hSigPureSA = new TH2F("hSigP_ITSPureSA",
 616                               "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
 617                               vX->GetNrows()-1,vX->GetMatrixArray(),
 618                               300,0,300);
 619   fListQAitsPureSA->Add(hSigPureSA);
 620
 621   delete vX;
 622 }
 623
 624 //______________________________________________________________________________
 625 void AliAnalysisTaskPIDqa::SetupTPCqa()
 626 {
 627   //
 628   // Create the TPC qa objects
 629   //
 630
 631   TVectorD *vX=MakeLogBinning(200,.1,30);
 632
 633   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 634     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_%s",AliPID::ParticleName(ispecie)),
 635                               Form("TPC n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 636                               vX->GetNrows()-1,vX->GetMatrixArray(),
 637                               200,-10,10);
 638     fListQAtpc->Add(hNsigmaP);
 639   }
 640
 641
 642   TH2F *hSig = new TH2F("hSigP_TPC",
 643                         "TPC signal vs. p;p [GeV]; TPC signal [arb. units]",
 644                         vX->GetNrows()-1,vX->GetMatrixArray(),
 645                         300,0,300);
 646   fListQAtpc->Add(hSig);
 647
 648   delete vX;
 649 }
 650
 651 //______________________________________________________________________________
 652 void AliAnalysisTaskPIDqa::SetupTRDqa()
 653 {
 654   //
 655   // Create the TRD qa objects
 656   //
 657   TVectorD *vX=MakeLogBinning(200,.1,30);
 658   for(Int_t itl = 0; itl < 6; ++itl){
 659     for(Int_t ispecie = 0; ispecie < AliPID::kSPECIES; ispecie++){
 660       TH2F *hLikeP = new TH2F(Form("hLikeP_TRD_%dtls_%s", itl, AliPID::ParticleName(ispecie)),
 661                               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)),
 662                               vX->GetNrows()-1, vX->GetMatrixArray(),
 663                               100, 0., 1.);
 664       fListQAtrd->Add(hLikeP);
 665     }
 666   }
 667   delete vX;
 668 }
 669
 670 //______________________________________________________________________________
 671 void AliAnalysisTaskPIDqa::SetupTOFqa()
 672 {
 673   //
 674   // Create the TOF qa objects
 675   //
 676
 677   TVectorD *vX=MakeLogBinning(200,.1,30);
 678
 679   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 680     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_%s",AliPID::ParticleName(ispecie)),
 681                               Form("TOF n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 682                               vX->GetNrows()-1,vX->GetMatrixArray(),
 683                               200,-10,10);
 684     fListQAtof->Add(hNsigmaP);
 685   }
 686
 687   // for Kaons PID we differentiate on Time Zero
 688   TH1F *hnSigT0Fill = new TH1F("hNsigma_TOF_Kaon_T0-Fill","TOF n#sigma (Kaon) T0-FILL [1-2. GeV/c]",200,-10,10);
 689   fListQAtof->Add(hnSigT0Fill);
 690   TH1F *hnSigT0T0 = new TH1F("hNsigma_TOF_Kaon_T0-T0","TOF n#sigma (Kaon) T0-T0 [1-2. GeV/c]",200,-10,10);
 691   fListQAtof->Add(hnSigT0T0);
 692   TH1F *hnSigT0TOF = new TH1F("hNsigma_TOF_Kaon_T0-TOF","TOF n#sigma (Kaon) T0-TOF [1.-2. GeV/c]",200,-10,10);
 693   fListQAtof->Add(hnSigT0TOF);
 694   TH1F *hnSigT0Best = new TH1F("hNsigma_TOF_Kaon_T0-Best","TOF n#sigma (Kaon) T0-Best [1-2. GeV/c]",200,-10,10);
 695   fListQAtof->Add(hnSigT0Best);
 696
 697   TH2F *hSig = new TH2F("hSigP_TOF",
 698                         "TOF signal vs. p;p [GeV]; TOF signal [arb. units]",
 699                         vX->GetNrows()-1,vX->GetMatrixArray(),
 700                         300,0,300);
 701
 702   delete vX;
 703
 704   fListQAtof->Add(hSig);
 705
 706   TH1F *hStartTimeMaskTOF = new TH1F("hStartTimeMask_TOF","StartTime mask",8,0,8);
 707   fListQAtof->Add(hStartTimeMaskTOF);
 708   TH1F *hStartTimeResTOF = new TH1F("hStartTimeRes_TOF","StartTime resolution [ps]",100,0,500);
 709   fListQAtof->Add(hStartTimeResTOF);
 710
 711   TH1F *hnTracksAtTOF = new TH1F("hnTracksAt_TOF","Matched tracks at TOF",20,0,20);
 712   fListQAtof->Add(hnTracksAtTOF);
 713   TH1F *hT0MakerEff = new TH1F("hT0MakerEff","Events with T0-TOF vs nTracks",20,0,20);
 714   fListQAtof->Add(hT0MakerEff);
 715
 716   // this in principle should stay on a T0 PID QA, but are just the data prepared for TOF use
 717   TH1F *hStartTimeAT0 = new TH1F("hStartTimeA_T0","StartTime from T0A [ps]",1000,-1000,1000);
 718   fListQAtof->Add(hStartTimeAT0);
 719   TH1F *hStartTimeCT0 = new TH1F("hStartTimeC_T0","StartTime from T0C [ps]",1000,-1000,1000);
 720   fListQAtof->Add(hStartTimeCT0);
 721   TH1F *hStartTimeACT0 = new TH1F("hStartTimeAC_T0","StartTime from T0AC [ps]",1000,-1000,1000);;
 722   fListQAtof->Add(hStartTimeACT0);
 723 }
 724
 725 //______________________________________________________________________________
 726 void AliAnalysisTaskPIDqa::SetupEMCALqa()
 727 {
 728   //
 729   // Create the EMCAL qa objects
 730   //
 731
 732   TVectorD *vX=MakeLogBinning(200,.1,30);
 733
 734   TH2F *hNsigmaPt = new TH2F(Form("hNsigmaPt_EMCAL_%s",AliPID::ParticleName(0)),
 735                              Form("EMCAL n#sigma %s vs. p_{T};p_{T} [GeV]; n#sigma",AliPID::ParticleName(0)),
 736                              vX->GetNrows()-1,vX->GetMatrixArray(),
 737                              200,-10,10);
 738   fListQAemcal->Add(hNsigmaPt);
 739
 740   TH2F *hSigPt = new TH2F("hSigPt_EMCAL",
 741                         "EMCAL signal (E/p) vs. p_{T};p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
 742                         vX->GetNrows()-1,vX->GetMatrixArray(),
 743                         200,0,2);
 744   fListQAemcal->Add(hSigPt);
 745
 746   delete vX;
 747 }
 748
 749 //______________________________________________________________________________
 750 void AliAnalysisTaskPIDqa::SetupTPCTOFqa()
 751 {
 752   //
 753   // Create the qa objects for TPC + TOF combination
 754   //
 755
 756   TVectorD *vX=MakeLogBinning(200,.1,30);
 757
 758   //TPC signals after TOF cut
 759   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 760     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_TOF_%s",AliPID::ParticleName(ispecie)),
 761                               Form("TPC n#sigma %s vs. p (after TOF 3#sigma cut);p_{TPC} [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 762                               vX->GetNrows()-1,vX->GetMatrixArray(),
 763                               200,-10,10);
 764     fListQAtpctof->Add(hNsigmaP);
 765   }
 766
 767   //TOF signals after TPC cut
 768   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 769     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_TPC_%s",AliPID::ParticleName(ispecie)),
 770                               Form("TOF n#sigma %s vs. p (after TPC n#sigma cut);p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 771                               vX->GetNrows()-1,vX->GetMatrixArray(),
 772                               200,-10,10);
 773     fListQAtpctof->Add(hNsigmaP);
 774   }
 775
 776   //EMCAL signal after TOF and TPC cut
 777   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 778     TH2F *heopPt = new TH2F(Form("heopPt_TOF_TPC_%s",AliPID::ParticleName(ispecie)),
 779                             Form("EMCAL signal (E/p) %s vs. p_{T};p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",AliPID::ParticleName(ispecie)),
 780                             vX->GetNrows()-1,vX->GetMatrixArray(),
 781                             200,0,2);
 782     fListQAtpctof->Add(heopPt);
 783   }
 784
 785   delete vX;
 786 }
 787
 788 //______________________________________________________________________________
 789 TVectorD* AliAnalysisTaskPIDqa::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
 790 {
 791   //
 792   // Make logarithmic binning
 793   // the user has to delete the array afterwards!!!
 794   //
 795
 796   //check limits
 797   if (xmin<1e-20 || xmax<1e-20){
 798     AliError("For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
 799     return MakeLinBinning(nbinsX, xmin, xmax);
 800   }
 801   if (xmax<xmin){
 802     Double_t tmp=xmin;
 803     xmin=xmax;
 804     xmax=tmp;
 805   }
 806   TVectorD *binLim=new TVectorD(nbinsX+1);
 807   Double_t first=xmin;
 808   Double_t last=xmax;
 809   Double_t expMax=TMath::Log(last/first);
 810   for (Int_t i=0; i<nbinsX+1; ++i){
 811     (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
 812   }
 813   return binLim;
 814 }
 815
 816 //______________________________________________________________________________
 817 TVectorD* AliAnalysisTaskPIDqa::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
 818 {
 819   //
 820   // Make linear binning
 821   // the user has to delete the array afterwards!!!
 822   //
 823   if (xmax<xmin){
 824     Double_t tmp=xmin;
 825     xmin=xmax;
 826     xmax=tmp;
 827   }
 828   TVectorD *binLim=new TVectorD(nbinsX+1);
 829   Double_t first=xmin;
 830   Double_t last=xmax;
 831   Double_t binWidth=(last-first)/nbinsX;
 832   for (Int_t i=0; i<nbinsX+1; ++i){
 833     (*binLim)[i]=first+binWidth*(Double_t)i;
 834   }
 835   return binLim;
 836 }
 837
 838 //_____________________________________________________________________________
 839 TVectorD* AliAnalysisTaskPIDqa::MakeArbitraryBinning(const char* bins)
 840 {
 841   //
 842   // Make arbitrary binning, bins separated by a ','
 843   //
 844   TString limits(bins);
 845   if (limits.IsNull()){
 846     AliError("Bin Limit string is empty, cannot add the variable");
 847     return 0x0;
 848   }
 849
 850   TObjArray *arr=limits.Tokenize(",");
 851   Int_t nLimits=arr->GetEntries();
 852   if (nLimits<2){
 853     AliError("Need at leas 2 bin limits, cannot add the variable");
 854     delete arr;
 855     return 0x0;
 856   }
 857
 858   TVectorD *binLimits=new TVectorD(nLimits);
 859   for (Int_t iLim=0; iLim<nLimits; ++iLim){
 860     (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
 861   }
 862
 863   delete arr;
 864   return binLimits;
 865 }
 866