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