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       else{
 482         Printf("status status = AliVTrack::kEMCALmatch, BUT no matched cluster!");
 483       }
 484     }
 485   }
 486 }
 487
 488 //______________________________________________________________________________
 489 void AliAnalysisTaskPIDqa::FillTPCTOFqa()
 490 {
 491   //
 492   // Fill PID qa histograms for the TOF
 493   //   Here also the TPC histograms after TOF selection are filled
 494   //
 495
 496   AliVEvent *event=InputEvent();
 497
 498   Int_t ntracks=event->GetNumberOfTracks();
 499   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 500     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 501
 502     //
 503     //basic track cuts
 504     //
 505     ULong_t status=track->GetStatus();
 506     // not that nice. status bits not in virtual interface
 507     // TPC refit + ITS refit +
 508     // TOF out + TOFpid +
 509     // kTIME
 510     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 511         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
 512         !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid ) ||
 513         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||
 514         !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||
 515         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;
 516
 517     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 518     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 519     if (track->GetTPCNclsF()>0) {
 520       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 521     }
 522
 523     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 524
 525
 526     Double_t mom=track->P();
 527     Double_t momTPC=track->GetTPCmomentum();
 528
 529     for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 530       //TOF nSigma
 531       Double_t nSigmaTOF=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
 532       Double_t nSigmaTPC=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
 533
 534       //TPC after TOF cut
 535       TH2 *h=(TH2*)fListQAtpctof->At(ispecie);
 536       if (h && TMath::Abs(nSigmaTOF)<3.) h->Fill(momTPC,nSigmaTPC);
 537
 538       //TOF after TPC cut
 539       h=(TH2*)fListQAtpctof->At(ispecie+AliPID::kSPECIES);
 540       if (h && TMath::Abs(nSigmaTPC)<3.) h->Fill(mom,nSigmaTOF);
 541
 542       //EMCAL after TOF and TPC cut
 543       h=(TH2*)fListQAtpctof->At(ispecie+2*AliPID::kSPECIES);
 544       if (h && TMath::Abs(nSigmaTOF)<3. && TMath::Abs(nSigmaTPC)<3. ){
 545
 546         Int_t nMatchClus = track->GetEMCALcluster();
 547         Double_t pt      = track->Pt();
 548         Double_t eop     = -1.;
 549
 550         if(nMatchClus > -1){
 551
 552           AliVCluster *matchedClus = (AliVCluster*)event->GetCaloCluster(nMatchClus);
 553
 554           if(matchedClus){
 555
 556             // matched cluster is EMCAL
 557             if(matchedClus->IsEMCAL()){
 558
 559               Double_t fClsE       = matchedClus->E();
 560               eop                  = fClsE/mom;
 561
 562               h->Fill(pt,eop);
 563
 564
 565             }
 566           }
 567         }
 568       }
 569     }
 570   }
 571 }
 572
 573 //______________________________________________________________________________
 574 void AliAnalysisTaskPIDqa::SetupITSqa()
 575 {
 576   //
 577   // Create the ITS qa objects
 578   //
 579
 580   TVectorD *vX=MakeLogBinning(200,.1,30);
 581
 582   //ITS+TPC tracks
 583   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 584     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_ITS_%s",AliPID::ParticleName(ispecie)),
 585                               Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 586                               vX->GetNrows()-1,vX->GetMatrixArray(),
 587                               200,-10,10);
 588     fListQAits->Add(hNsigmaP);
 589   }
 590   TH2F *hSig = new TH2F("hSigP_ITS",
 591                         "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
 592                         vX->GetNrows()-1,vX->GetMatrixArray(),
 593                         300,0,300);
 594   fListQAits->Add(hSig);
 595
 596   //ITS Standalone tracks
 597   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 598     TH2F *hNsigmaPSA = new TH2F(Form("hNsigmaP_ITSSA_%s",AliPID::ParticleName(ispecie)),
 599                                 Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 600                                 vX->GetNrows()-1,vX->GetMatrixArray(),
 601                                 200,-10,10);
 602     fListQAitsSA->Add(hNsigmaPSA);
 603   }
 604   TH2F *hSigSA = new TH2F("hSigP_ITSSA",
 605                           "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
 606                           vX->GetNrows()-1,vX->GetMatrixArray(),
 607                           300,0,300);
 608   fListQAitsSA->Add(hSigSA);
 609
 610   //ITS Pure Standalone tracks
 611   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 612     TH2F *hNsigmaPPureSA = new TH2F(Form("hNsigmaP_ITSPureSA_%s",AliPID::ParticleName(ispecie)),
 613                                     Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 614                                     vX->GetNrows()-1,vX->GetMatrixArray(),
 615                                     200,-10,10);
 616     fListQAitsPureSA->Add(hNsigmaPPureSA);
 617   }
 618   TH2F *hSigPureSA = new TH2F("hSigP_ITSPureSA",
 619                               "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
 620                               vX->GetNrows()-1,vX->GetMatrixArray(),
 621                               300,0,300);
 622   fListQAitsPureSA->Add(hSigPureSA);
 623
 624   delete vX;
 625 }
 626
 627 //______________________________________________________________________________
 628 void AliAnalysisTaskPIDqa::SetupTPCqa()
 629 {
 630   //
 631   // Create the TPC qa objects
 632   //
 633
 634   TVectorD *vX=MakeLogBinning(200,.1,30);
 635
 636   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 637     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_%s",AliPID::ParticleName(ispecie)),
 638                               Form("TPC n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 639                               vX->GetNrows()-1,vX->GetMatrixArray(),
 640                               200,-10,10);
 641     fListQAtpc->Add(hNsigmaP);
 642   }
 643
 644
 645   TH2F *hSig = new TH2F("hSigP_TPC",
 646                         "TPC signal vs. p;p [GeV]; TPC signal [arb. units]",
 647                         vX->GetNrows()-1,vX->GetMatrixArray(),
 648                         300,0,300);
 649   fListQAtpc->Add(hSig);
 650
 651   delete vX;
 652 }
 653
 654 //______________________________________________________________________________
 655 void AliAnalysisTaskPIDqa::SetupTRDqa()
 656 {
 657   //
 658   // Create the TRD qa objects
 659   //
 660   TVectorD *vX=MakeLogBinning(200,.1,30);
 661   for(Int_t itl = 0; itl < 6; ++itl){
 662     for(Int_t ispecie = 0; ispecie < AliPID::kSPECIES; ispecie++){
 663       TH2F *hLikeP = new TH2F(Form("hLikeP_TRD_%dtls_%s", itl, AliPID::ParticleName(ispecie)),
 664                               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)),
 665                               vX->GetNrows()-1, vX->GetMatrixArray(),
 666                               100, 0., 1.);
 667       fListQAtrd->Add(hLikeP);
 668     }
 669   }
 670   delete vX;
 671 }
 672
 673 //______________________________________________________________________________
 674 void AliAnalysisTaskPIDqa::SetupTOFqa()
 675 {
 676   //
 677   // Create the TOF qa objects
 678   //
 679
 680   TVectorD *vX=MakeLogBinning(200,.1,30);
 681
 682   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 683     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_%s",AliPID::ParticleName(ispecie)),
 684                               Form("TOF n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 685                               vX->GetNrows()-1,vX->GetMatrixArray(),
 686                               200,-10,10);
 687     fListQAtof->Add(hNsigmaP);
 688   }
 689
 690   // for Kaons PID we differentiate on Time Zero
 691   TH1F *hnSigT0Fill = new TH1F("hNsigma_TOF_Kaon_T0-Fill","TOF n#sigma (Kaon) T0-FILL [1-2. GeV/c]",200,-10,10);
 692   fListQAtof->Add(hnSigT0Fill);
 693   TH1F *hnSigT0T0 = new TH1F("hNsigma_TOF_Kaon_T0-T0","TOF n#sigma (Kaon) T0-T0 [1-2. GeV/c]",200,-10,10);
 694   fListQAtof->Add(hnSigT0T0);
 695   TH1F *hnSigT0TOF = new TH1F("hNsigma_TOF_Kaon_T0-TOF","TOF n#sigma (Kaon) T0-TOF [1.-2. GeV/c]",200,-10,10);
 696   fListQAtof->Add(hnSigT0TOF);
 697   TH1F *hnSigT0Best = new TH1F("hNsigma_TOF_Kaon_T0-Best","TOF n#sigma (Kaon) T0-Best [1-2. GeV/c]",200,-10,10);
 698   fListQAtof->Add(hnSigT0Best);
 699
 700
 701   TH2F *hSig = new TH2F("hSigP_TOF",
 702                         "TOF signal vs. p;p [GeV]; TOF signal [arb. units]",
 703                         vX->GetNrows()-1,vX->GetMatrixArray(),
 704                         300,0,300);
 705
 706   delete vX;
 707
 708   fListQAtof->Add(hSig);
 709
 710   TH1F *hStartTimeMaskTOF = new TH1F("hStartTimeMask_TOF","StartTime mask",8,0,8);
 711   fListQAtof->Add(hStartTimeMaskTOF);
 712   TH1F *hStartTimeResTOF = new TH1F("hStartTimeRes_TOF","StartTime resolution [ps]",100,0,500);
 713   fListQAtof->Add(hStartTimeResTOF);
 714
 715   TH1F *hnTracksAtTOF = new TH1F("hnTracksAt_TOF","Matched tracks at TOF",20,0,20);
 716   fListQAtof->Add(hnTracksAtTOF);
 717   TH1F *hT0MakerEff = new TH1F("hT0MakerEff","Events with T0-TOF vs nTracks",20,0,20);
 718   fListQAtof->Add(hT0MakerEff);
 719
 720   // this in principle should stay on a T0 PID QA, but are just the data prepared for TOF use
 721   TH1F *hStartTimeAT0 = new TH1F("hStartTimeA_T0","StartTime from T0A [ps]",1000,-1000,1000);
 722   fListQAtof->Add(hStartTimeAT0);
 723   TH1F *hStartTimeCT0 = new TH1F("hStartTimeC_T0","StartTime from T0C [ps]",1000,-1000,1000);
 724   fListQAtof->Add(hStartTimeCT0);
 725   TH1F *hStartTimeACT0 = new TH1F("hStartTimeAC_T0","StartTime from T0AC [ps]",1000,-1000,1000);;
 726   fListQAtof->Add(hStartTimeACT0);
 727 }
 728
 729 //______________________________________________________________________________
 730 void AliAnalysisTaskPIDqa::SetupEMCALqa()
 731 {
 732   //
 733   // Create the EMCAL qa objects
 734   //
 735
 736   TVectorD *vX=MakeLogBinning(200,.1,30);
 737
 738   TH2F *hNsigmaPt = new TH2F(Form("hNsigmaPt_EMCAL_%s",AliPID::ParticleName(0)),
 739                              Form("EMCAL n#sigma %s vs. p_{T};p_{T} [GeV]; n#sigma",AliPID::ParticleName(0)),
 740                              vX->GetNrows()-1,vX->GetMatrixArray(),
 741                              200,-10,10);
 742   fListQAemcal->Add(hNsigmaPt);
 743
 744   TH2F *hSigPt = new TH2F("hSigPt_EMCAL",
 745                         "EMCAL signal (E/p) vs. p_{T};p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",
 746                         vX->GetNrows()-1,vX->GetMatrixArray(),
 747                         200,0,2);
 748   fListQAemcal->Add(hSigPt);
 749
 750   delete vX;
 751 }
 752
 753 //______________________________________________________________________________
 754 void AliAnalysisTaskPIDqa::SetupTPCTOFqa()
 755 {
 756   //
 757   // Create the qa objects for TPC + TOF combination
 758   //
 759
 760   TVectorD *vX=MakeLogBinning(200,.1,30);
 761
 762   //TPC signals after TOF cut
 763   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 764     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_TOF_%s",AliPID::ParticleName(ispecie)),
 765                               Form("TPC n#sigma %s vs. p (after TOF 3#sigma cut);p_{TPC} [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 766                               vX->GetNrows()-1,vX->GetMatrixArray(),
 767                               200,-10,10);
 768     fListQAtpctof->Add(hNsigmaP);
 769   }
 770
 771   //TOF signals after TPC cut
 772   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 773     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_TPC_%s",AliPID::ParticleName(ispecie)),
 774                               Form("TOF n#sigma %s vs. p (after TPC n#sigma cut);p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 775                               vX->GetNrows()-1,vX->GetMatrixArray(),
 776                               200,-10,10);
 777     fListQAtpctof->Add(hNsigmaP);
 778   }
 779
 780   //EMCAL signal after TOF and TPC cut
 781   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 782     TH2F *heopPt = new TH2F(Form("heopPt_TOF_TPC_%s",AliPID::ParticleName(ispecie)),
 783                             Form("EMCAL signal (E/p) %s vs. p_{T};p_{T} [GeV]; EMCAL signal (E/p) [arb. units]",AliPID::ParticleName(ispecie)),
 784                             vX->GetNrows()-1,vX->GetMatrixArray(),
 785                             200,0,2);
 786     fListQAtpctof->Add(heopPt);
 787   }
 788
 789   delete vX;
 790 }
 791
 792 //______________________________________________________________________________
 793 TVectorD* AliAnalysisTaskPIDqa::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
 794 {
 795   //
 796   // Make logarithmic binning
 797   // the user has to delete the array afterwards!!!
 798   //
 799
 800   //check limits
 801   if (xmin<1e-20 || xmax<1e-20){
 802     AliError("For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
 803     return MakeLinBinning(nbinsX, xmin, xmax);
 804   }
 805   if (xmax<xmin){
 806     Double_t tmp=xmin;
 807     xmin=xmax;
 808     xmax=tmp;
 809   }
 810   TVectorD *binLim=new TVectorD(nbinsX+1);
 811   Double_t first=xmin;
 812   Double_t last=xmax;
 813   Double_t expMax=TMath::Log(last/first);
 814   for (Int_t i=0; i<nbinsX+1; ++i){
 815     (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
 816   }
 817   return binLim;
 818 }
 819
 820 //______________________________________________________________________________
 821 TVectorD* AliAnalysisTaskPIDqa::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
 822 {
 823   //
 824   // Make linear binning
 825   // the user has to delete the array afterwards!!!
 826   //
 827   if (xmax<xmin){
 828     Double_t tmp=xmin;
 829     xmin=xmax;
 830     xmax=tmp;
 831   }
 832   TVectorD *binLim=new TVectorD(nbinsX+1);
 833   Double_t first=xmin;
 834   Double_t last=xmax;
 835   Double_t binWidth=(last-first)/nbinsX;
 836   for (Int_t i=0; i<nbinsX+1; ++i){
 837     (*binLim)[i]=first+binWidth*(Double_t)i;
 838   }
 839   return binLim;
 840 }
 841
 842 //_____________________________________________________________________________
 843 TVectorD* AliAnalysisTaskPIDqa::MakeArbitraryBinning(const char* bins)
 844 {
 845   //
 846   // Make arbitrary binning, bins separated by a ','
 847   //
 848   TString limits(bins);
 849   if (limits.IsNull()){
 850     AliError("Bin Limit string is empty, cannot add the variable");
 851     return 0x0;
 852   }
 853
 854   TObjArray *arr=limits.Tokenize(",");
 855   Int_t nLimits=arr->GetEntries();
 856   if (nLimits<2){
 857     AliError("Need at leas 2 bin limits, cannot add the variable");
 858     delete arr;
 859     return 0x0;
 860   }
 861
 862   TVectorD *binLimits=new TVectorD(nLimits);
 863   for (Int_t iLim=0; iLim<nLimits; ++iLim){
 864     (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
 865   }
 866
 867   delete arr;
 868   return binLimits;
 869 }
 870