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   //
 338   // Fill PID qa histograms for the TOF
 339   //
 340
 341   AliVEvent *event=InputEvent();
 342
 343
 344   Int_t ntracks=event->GetNumberOfTracks();
 345   Int_t tracksAtTof = 0;
 346   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 347     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 348
 349     //
 350     //basic track cuts
 351     //
 352     ULong_t status=track->GetStatus();
 353     // TPC refit + ITS refit +
 354     // TOF out + TOFpid +
 355     // kTIME
 356     // (we don't use kTOFmismatch because it depends on TPC....)
 357     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 358         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
 359         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||
 360         !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||
 361         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;
 362
 363     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 364     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 365     if (track->GetTPCNclsF()>0) {
 366       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 367     }
 368
 369     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 370
 371     tracksAtTof++;
 372
 373     Double_t mom=track->P();
 374
 375     for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 376       TH2 *h=(TH2*)fListQAtof->At(ispecie);
 377       if (!h) continue;
 378       Double_t nSigma=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
 379       h->Fill(mom,nSigma);
 380     }
 381
 382     TH2 *h=(TH2*)fListQAtof->FindObject("hSigP_TOF");
 383     if (h) {
 384       Double_t sig=track->GetTOFsignal();
 385       h->Fill(mom,sig);
 386     }
 387
 388     Double_t mask = (Double_t)fPIDResponse->GetTOFResponse().GetStartTimeMask(mom) + 0.5;
 389     ((TH1F*)fListQAtof->FindObject("hStartTimeMask_TOF"))->Fill(mask);
 390
 391     Double_t res = (Double_t)fPIDResponse->GetTOFResponse().GetStartTimeRes(mom);
 392     ((TH1F*)fListQAtof->FindObject("hStartTimeRes_TOF"))->Fill(res);
 393
 394     AliESDEvent *esd = dynamic_cast<AliESDEvent *>(event);
 395     if (esd) {
 396       Double_t startTime = esd->GetT0TOF(0);
 397       if (startTime < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeAC_T0"))->Fill(startTime);
 398       else {
 399         startTime = esd->GetT0TOF(1);
 400         if (startTime < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeA_T0"))->Fill(startTime);
 401         startTime = esd->GetT0TOF(2);
 402         if (startTime < 90000) ((TH1F*)fListQAtof->FindObject("hStartTimeC_T0"))->Fill(startTime);
 403       }
 404     }
 405   }
 406   if (tracksAtTof > 0) {
 407     ((TH1F* )fListQAtof->FindObject("hnTracksAt_TOF"))->Fill(tracksAtTof);
 408     Int_t mask = fPIDResponse->GetTOFResponse().GetStartTimeMask(5.);
 409     if (mask & 0x1) ((TH1F*)fListQAtof->FindObject("hT0MakerEff"))->Fill(tracksAtTof);
 410   }
 411
 412 }
 413
 414 //______________________________________________________________________________
 415 void AliAnalysisTaskPIDqa::FillEMCALqa()
 416 {
 417   //
 418   // Fill PID qa histograms for the EMCAL
 419   //
 420
 421   AliVEvent *event=InputEvent();
 422
 423   Int_t ntracks=event->GetNumberOfTracks();
 424   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 425     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 426
 427     //
 428     //basic track cuts
 429     //
 430     ULong_t status=track->GetStatus();
 431     // not that nice. status bits not in virtual interface
 432     // TPC refit + ITS refit +
 433     // TOF out + TOFpid +
 434     // kTIME
 435     if (!( (status & AliVTrack::kEMCALmatch) == AliVTrack::kEMCALmatch) ) continue;
 436
 437     Double_t pt=track->Pt();
 438
 439     //EMCAL nSigma (only for electrons at the moment)
 440     TH2 *h=(TH2*)fListQAemcal->At(0);
 441     if (!h) continue;
 442     Double_t nSigma=fPIDResponse->NumberOfSigmasEMCAL(track, (AliPID::EParticleType)0);
 443     h->Fill(pt,nSigma);
 444
 445     //EMCAL signal (E/p vs. pT)
 446     h=(TH2*)fListQAemcal->At(1);
 447     if (h) {
 448
 449       Int_t nMatchClus = track->GetEMCALcluster();
 450       Double_t mom     = track->P();
 451       Double_t eop     = -1.;
 452
 453       if(nMatchClus > -1){
 454
 455         AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
 456
 457         if(matchedClus){
 458
 459           // matched cluster is EMCAL
 460           if(matchedClus->IsEMCAL()){
 461
 462             Double_t fClsE       = matchedClus->E();
 463             eop                  = fClsE/mom;
 464
 465             h->Fill(pt,eop);
 466
 467           }
 468         }
 469       }
 470       else{
 471         Printf("status status = AliVTrack::kEMCALmatch, BUT no matched cluster!");
 472       }
 473     }
 474   }
 475 }
 476
 477 //______________________________________________________________________________
 478 void AliAnalysisTaskPIDqa::FillTPCTOFqa()
 479 {
 480   //
 481   // Fill PID qa histograms for the TOF
 482   //   Here also the TPC histograms after TOF selection are filled
 483   //
 484
 485   AliVEvent *event=InputEvent();
 486
 487   Int_t ntracks=event->GetNumberOfTracks();
 488   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 489     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 490
 491     //
 492     //basic track cuts
 493     //
 494     ULong_t status=track->GetStatus();
 495     // not that nice. status bits not in virtual interface
 496     // TPC refit + ITS refit +
 497     // TOF out + TOFpid +
 498     // kTIME
 499     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 500         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
 501         !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid ) ||
 502         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||
 503         !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||
 504         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;
 505
 506     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 507     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 508     if (track->GetTPCNclsF()>0) {
 509       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 510     }
 511
 512     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 513
 514
 515     Double_t mom=track->P();
 516     Double_t momTPC=track->GetTPCmomentum();
 517
 518     for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 519       //TOF nSigma
 520       Double_t nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
 521       Double_t nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
 522
 523       //TPC after TOF cut
 524       TH2 *h=(TH2*)fListQAtpctof->At(ispecie);
 525       if (h && TMath::Abs(nSigmaTOF)<3.) h->Fill(momTPC,nSigmaTPC);
 526
 527       //TOF after TPC cut
 528       h=(TH2*)fListQAtpctof->At(ispecie+AliPID::kSPECIES);
 529       if (h && TMath::Abs(nSigmaTPC)<3.) h->Fill(mom,nSigmaTOF);
 530
 531       //EMCAL after TOF and TPC cut
 532       h=(TH2*)fListQAtpctof->At(ispecie+2*AliPID::kSPECIES);
 533       if (h && TMath::Abs(nSigmaTOF)<3. && TMath::Abs(nSigmaTPC)<3. ){
 534
 535         Int_t nMatchClus = track->GetEMCALcluster();
 536         Double_t pt      = track->Pt();
 537         Double_t eop     = -1.;
 538
 539         if(nMatchClus > -1){
 540
 541           AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
 542
 543           if(matchedClus){
 544
 545             // matched cluster is EMCAL
 546             if(matchedClus->IsEMCAL()){
 547
 548               Double_t fClsE       = matchedClus->E();
 549               eop                  = fClsE/mom;
 550
 551               h->Fill(pt,eop);
 552
 553
 554             }
 555           }
 556         }
 557       }
 558     }
 559   }
 560 }
 561
 562 //______________________________________________________________________________
 563 void AliAnalysisTaskPIDqa::SetupITSqa()
 564 {
 565   //
 566   // Create the ITS qa objects
 567   //
 568
 569   TVectorD *vX=MakeLogBinning(200,.1,30);
 570
 571   //ITS+TPC tracks
 572   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 573     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_ITS_%s",AliPID::ParticleName(ispecie)),
 574                               Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 575                               vX->GetNrows()-1,vX->GetMatrixArray(),
 576                               200,-10,10);
 577     fListQAits->Add(hNsigmaP);
 578   }
 579   TH2F *hSig = new TH2F("hSigP_ITS",
 580                         "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
 581                         vX->GetNrows()-1,vX->GetMatrixArray(),
 582                         300,0,300);
 583   fListQAits->Add(hSig);
 584
 585   //ITS Standalone tracks
 586   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 587     TH2F *hNsigmaPSA = new TH2F(Form("hNsigmaP_ITSSA_%s",AliPID::ParticleName(ispecie)),
 588                                 Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 589                                 vX->GetNrows()-1,vX->GetMatrixArray(),
 590                                 200,-10,10);
 591     fListQAitsSA->Add(hNsigmaPSA);
 592   }
 593   TH2F *hSigSA = new TH2F("hSigP_ITSSA",
 594                           "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
 595                           vX->GetNrows()-1,vX->GetMatrixArray(),
 596                           300,0,300);
 597   fListQAitsSA->Add(hSigSA);
 598
 599   //ITS Pure Standalone tracks
 600   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 601     TH2F *hNsigmaPPureSA = new TH2F(Form("hNsigmaP_ITSPureSA_%s",AliPID::ParticleName(ispecie)),
 602                                     Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 603                                     vX->GetNrows()-1,vX->GetMatrixArray(),
 604                                     200,-10,10);
 605     fListQAitsPureSA->Add(hNsigmaPPureSA);
 606   }
 607   TH2F *hSigPureSA = new TH2F("hSigP_ITSPureSA",
 608                               "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
 609                               vX->GetNrows()-1,vX->GetMatrixArray(),
 610                               300,0,300);
 611   fListQAitsPureSA->Add(hSigPureSA);
 612
 613   delete vX;
 614 }
 615
 616 //______________________________________________________________________________
 617 void AliAnalysisTaskPIDqa::SetupTPCqa()
 618 {
 619   //
 620   // Create the TPC qa objects
 621   //
 622
 623   TVectorD *vX=MakeLogBinning(200,.1,30);
 624
 625   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 626     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_%s",AliPID::ParticleName(ispecie)),
 627                               Form("TPC n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 628                               vX->GetNrows()-1,vX->GetMatrixArray(),
 629                               200,-10,10);
 630     fListQAtpc->Add(hNsigmaP);
 631   }
 632
 633
 634   TH2F *hSig = new TH2F("hSigP_TPC",
 635                         "TPC signal vs. p;p [GeV]; TPC signal [arb. units]",
 636                         vX->GetNrows()-1,vX->GetMatrixArray(),
 637                         300,0,300);
 638   fListQAtpc->Add(hSig);
 639
 640   delete vX;
 641 }
 642
 643 //______________________________________________________________________________
 644 void AliAnalysisTaskPIDqa::SetupTRDqa()
 645 {
 646   //
 647   // Create the TRD qa objects
 648   //
 649   TVectorD *vX=MakeLogBinning(200,.1,30);
 650   for(Int_t itl = 0; itl < 6; ++itl){
 651     for(Int_t ispecie = 0; ispecie < AliPID::kSPECIES; ispecie++){
 652       TH2F *hLikeP = new TH2F(Form("hLikeP_TRD_%dtls_%s", itl, AliPID::ParticleName(ispecie)),
 653                               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)),
 654                               vX->GetNrows()-1, vX->GetMatrixArray(),
 655                               100, 0., 1.);
 656       fListQAtrd->Add(hLikeP);
 657     }
 658   }
 659   delete vX;
 660 }
 661
 662 //______________________________________________________________________________
 663 void AliAnalysisTaskPIDqa::SetupTOFqa()
 664 {
 665   //
 666   // Create the TOF qa objects
 667   //
 668
 669   TVectorD *vX=MakeLogBinning(200,.1,30);
 670
 671   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 672     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_%s",AliPID::ParticleName(ispecie)),
 673                               Form("TOF n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 674                               vX->GetNrows()-1,vX->GetMatrixArray(),
 675                               200,-10,10);
 676     fListQAtof->Add(hNsigmaP);
 677     // to be added: t-texp without StartTime subtraction
 678   }
 679
 680
 681   TH2F *hSig = new TH2F("hSigP_TOF",
 682                         "TOF signal vs. p;p [GeV]; TOF signal [arb. units]",
 683                         vX->GetNrows()-1,vX->GetMatrixArray(),
 684                         300,0,300);
 685
 686   delete vX;
 687
 688   fListQAtof->Add(hSig);
 689
 690   TH1F *hStartTimeMaskTOF = new TH1F("hStartTimeMask_TOF","StartTime mask",8,0,8);
 691   fListQAtof->Add(hStartTimeMaskTOF);
 692   TH1F *hStartTimeResTOF = new TH1F("hStartTimeRes_TOF","StartTime resolution [ps]",100,0,500);
 693   fListQAtof->Add(hStartTimeResTOF);
 694
 695   TH1F *hnTracksAtTOF = new TH1F("hnTracksAt_TOF","Matched tracks at TOF",20,0,20);
 696   fListQAtof->Add(hnTracksAtTOF);
 697   TH1F *hT0MakerEff = new TH1F("hT0MakerEff","Events with T0-TOF vs nTracks",20,0,20);
 698   fListQAtof->Add(hT0MakerEff);
 699
 700   // this in principle should stay on a T0 PID QA, but are just the data prepared for TOF use
 701   TH1F *hStartTimeAT0 = new TH1F("hStartTimeA_T0","StartTime from T0A [ps]",1000,-1000,1000);
 702   fListQAtof->Add(hStartTimeAT0);
 703   TH1F *hStartTimeCT0 = new TH1F("hStartTimeC_T0","StartTime from T0C [ps]",1000,-1000,1000);
 704   fListQAtof->Add(hStartTimeCT0);
 705   TH1F *hStartTimeACT0 = new TH1F("hStartTimeAC_T0","StartTime from T0AC [ps]",1000,-1000,1000);;
 706   fListQAtof->Add(hStartTimeACT0);
 707 }
 708
 709 //______________________________________________________________________________
 710 void AliAnalysisTaskPIDqa::SetupEMCALqa()
 711 {
 712   //
 713   // Create the EMCAL qa objects
 714   //
 715
 716   TVectorD *vX=MakeLogBinning(200,.1,30);
 717
 718   TH2F *hNsigmaPt = new TH2F(Form("hNsigmaPt_EMCAL_%s",AliPID::ParticleName(0)),
 719                              Form("EMCAL n#sigma %s vs. p_{T};p_{T} [GeV]; n#sigma",AliPID::ParticleName(0)),
 720                              vX->GetNrows()-1,vX->GetMatrixArray(),
 721                              200,-10,10);
 722   fListQAemcal->Add(hNsigmaPt);
 723
 724   TH2F *hSigPt = new TH2F("hSigPt_EMCAL",
 725                         "EMCAL signal (E/p) vs. p_{T};p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
 726                         vX->GetNrows()-1,vX->GetMatrixArray(),
 727                         200,0,2);
 728   fListQAemcal->Add(hSigPt);
 729
 730   delete vX;
 731 }
 732
 733 //______________________________________________________________________________
 734 void AliAnalysisTaskPIDqa::SetupTPCTOFqa()
 735 {
 736   //
 737   // Create the qa objects for TPC + TOF combination
 738   //
 739
 740   TVectorD *vX=MakeLogBinning(200,.1,30);
 741
 742   //TPC signals after TOF cut
 743   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 744     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_TOF_%s",AliPID::ParticleName(ispecie)),
 745                               Form("TPC n#sigma %s vs. p (after TOF 3#sigma cut);p_{TPC} [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 746                               vX->GetNrows()-1,vX->GetMatrixArray(),
 747                               200,-10,10);
 748     fListQAtpctof->Add(hNsigmaP);
 749   }
 750
 751   //TOF signals after TPC cut
 752   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 753     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_TPC_%s",AliPID::ParticleName(ispecie)),
 754                               Form("TOF n#sigma %s vs. p (after TPC n#sigma cut);p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 755                               vX->GetNrows()-1,vX->GetMatrixArray(),
 756                               200,-10,10);
 757     fListQAtpctof->Add(hNsigmaP);
 758   }
 759
 760   //EMCAL signal after TOF and TPC cut
 761   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 762     TH2F *heopPt = new TH2F(Form("heopPt_TOF_TPC_%s",AliPID::ParticleName(ispecie)),
 763                             Form("EMCAL signal (E/p) %s vs. p_{T};p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",AliPID::ParticleName(ispecie)),
 764                             vX->GetNrows()-1,vX->GetMatrixArray(),
 765                             200,0,2);
 766     fListQAtpctof->Add(heopPt);
 767   }
 768
 769   delete vX;
 770 }
 771
 772 //______________________________________________________________________________
 773 TVectorD* AliAnalysisTaskPIDqa::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
 774 {
 775   //
 776   // Make logarithmic binning
 777   // the user has to delete the array afterwards!!!
 778   //
 779
 780   //check limits
 781   if (xmin<1e-20 || xmax<1e-20){
 782     AliError("For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
 783     return MakeLinBinning(nbinsX, xmin, xmax);
 784   }
 785   if (xmax<xmin){
 786     Double_t tmp=xmin;
 787     xmin=xmax;
 788     xmax=tmp;
 789   }
 790   TVectorD *binLim=new TVectorD(nbinsX+1);
 791   Double_t first=xmin;
 792   Double_t last=xmax;
 793   Double_t expMax=TMath::Log(last/first);
 794   for (Int_t i=0; i<nbinsX+1; ++i){
 795     (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
 796   }
 797   return binLim;
 798 }
 799
 800 //______________________________________________________________________________
 801 TVectorD* AliAnalysisTaskPIDqa::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
 802 {
 803   //
 804   // Make linear binning
 805   // the user has to delete the array afterwards!!!
 806   //
 807   if (xmax<xmin){
 808     Double_t tmp=xmin;
 809     xmin=xmax;
 810     xmax=tmp;
 811   }
 812   TVectorD *binLim=new TVectorD(nbinsX+1);
 813   Double_t first=xmin;
 814   Double_t last=xmax;
 815   Double_t binWidth=(last-first)/nbinsX;
 816   for (Int_t i=0; i<nbinsX+1; ++i){
 817     (*binLim)[i]=first+binWidth*(Double_t)i;
 818   }
 819   return binLim;
 820 }
 821
 822 //_____________________________________________________________________________
 823 TVectorD* AliAnalysisTaskPIDqa::MakeArbitraryBinning(const char* bins)
 824 {
 825   //
 826   // Make arbitrary binning, bins separated by a ','
 827   //
 828   TString limits(bins);
 829   if (limits.IsNull()){
 830     AliError("Bin Limit string is empty, cannot add the variable");
 831     return 0x0;
 832   }
 833
 834   TObjArray *arr=limits.Tokenize(",");
 835   Int_t nLimits=arr->GetEntries();
 836   if (nLimits<2){
 837     AliError("Need at leas 2 bin limits, cannot add the variable");
 838     delete arr;
 839     return 0x0;
 840   }
 841
 842   TVectorD *binLimits=new TVectorD(nLimits);
 843   for (Int_t iLim=0; iLim<nLimits; ++iLim){
 844     (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
 845   }
 846
 847   delete arr;
 848   return binLimits;
 849 }
 850