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