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 "AliAnalysisTaskPIDqa.h"
  42
  43
  44 ClassImp(AliAnalysisTaskPIDqa)
  45
  46 //______________________________________________________________________________
  47 AliAnalysisTaskPIDqa::AliAnalysisTaskPIDqa():
  48 AliAnalysisTaskSE(),
  49 fPIDResponse(0x0),
  50 fListQA(0x0),
  51 fListQAits(0x0),
  52 fListQAtpc(0x0),
  53 fListQAtrd(0x0),
  54 fListQAtof(0x0)
  55 {
  56   //
  57   // Dummy constructor
  58   //
  59 }
  60
  61 //______________________________________________________________________________
  62 AliAnalysisTaskPIDqa::AliAnalysisTaskPIDqa(const char* name):
  63 AliAnalysisTaskSE(name),
  64 fPIDResponse(0x0),
  65 fListQA(0x0),
  66 fListQAits(0x0),
  67 fListQAtpc(0x0),
  68 fListQAtrd(0x0),
  69 fListQAtof(0x0)
  70 {
  71   //
  72   // Default constructor
  73   //
  74   DefineInput(0,TChain::Class());
  75   DefineOutput(1,TList::Class());
  76 }
  77
  78 //______________________________________________________________________________
  79 AliAnalysisTaskPIDqa::~AliAnalysisTaskPIDqa()
  80 {
  81   //
  82   // Destructor
  83   //
  84
  85 }
  86
  87 //______________________________________________________________________________
  88 void AliAnalysisTaskPIDqa::UserCreateOutputObjects()
  89 {
  90   //
  91   // Create the output QA objects
  92   //
  93
  94   AliLog::SetClassDebugLevel("AliAnalysisTaskPIDqa",10);
  95
  96   //input hander
  97   AliAnalysisManager *man=AliAnalysisManager::GetAnalysisManager();
  98   AliInputEventHandler *inputHandler=dynamic_cast<AliInputEventHandler*>(man->GetInputEventHandler());
  99   if (!inputHandler) AliFatal("Input handler needed");
 100
 101   //pid response object
 102   fPIDResponse=inputHandler->GetPIDResponse();
 103   if (!fPIDResponse) AliFatal("PIDResponse object was not created");
 104
 105   //
 106   fListQA=new TList;
 107   fListQA->SetOwner();
 108
 109   fListQAits=new TList;
 110   fListQAits->SetOwner();
 111   fListQAits->SetName("ITS");
 112
 113   fListQAtpc=new TList;
 114   fListQAtpc->SetOwner();
 115   fListQAtpc->SetName("TPC");
 116
 117   fListQAtrd=new TList;
 118   fListQAtrd->SetOwner();
 119   fListQAtrd->SetName("TRD");
 120
 121   fListQAtof=new TList;
 122   fListQAtof->SetOwner();
 123   fListQAtof->SetName("TOF");
 124
 125   fListQA->Add(fListQAits);
 126   fListQA->Add(fListQAtpc);
 127   fListQA->Add(fListQAtrd);
 128   fListQA->Add(fListQAtof);
 129
 130   SetupITSqa();
 131   SetupTPCqa();
 132   SetupTRDqa();
 133   SetupTOFqa();
 134
 135   PostData(1,fListQA);
 136 }
 137
 138
 139 //______________________________________________________________________________
 140 void AliAnalysisTaskPIDqa::UserExec(Option_t */*option*/)
 141 {
 142   //
 143   // Setup the PID response functions and fill the QA histograms
 144   //
 145
 146   AliVEvent *event=InputEvent();
 147   if (!event||!fPIDResponse) return;
 148
 149
 150   FillITSqa();
 151   FillTPCqa();
 152   FillTOFqa();
 153
 154   PostData(1,fListQA);
 155 }
 156
 157 //______________________________________________________________________________
 158 void AliAnalysisTaskPIDqa::FillITSqa()
 159 {
 160   AliVEvent *event=InputEvent();
 161
 162   Int_t ntracks=event->GetNumberOfTracks();
 163   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 164     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 165     ULong_t status=track->GetStatus();
 166     // not that nice. status bits not in virtual interface
 167     // ITS refit + ITS pid
 168     if (!( ( (status&0x0004)==0x0004 ) && ( (status&0x0008)==0x0008 ) )) return;
 169     Double_t mom=track->P();
 170
 171     for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 172       TH2 *h=(TH2*)fListQAits->At(ispecie);
 173       if (!h) continue;
 174       Double_t nSigma=fPIDResponse->NumberOfSigmasITS(track, (AliPID::EParticleType)ispecie);
 175       h->Fill(mom,nSigma);
 176     }
 177
 178     TH2 *h=(TH2*)fListQAits->At(AliPID::kSPECIES);
 179     if (h) {
 180       Double_t sig=track->GetITSsignal();
 181       h->Fill(mom,sig);
 182     }
 183
 184   }
 185 }
 186
 187 //______________________________________________________________________________
 188 void AliAnalysisTaskPIDqa::FillTPCqa()
 189 {
 190   AliVEvent *event=InputEvent();
 191
 192   Int_t ntracks=event->GetNumberOfTracks();
 193   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 194     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 195
 196     //
 197     //basic track cuts
 198     //
 199     ULong_t status=track->GetStatus();
 200     // not that nice. status bits not in virtual interface
 201     // TPC refit + ITS refit + TPC pid
 202     if (!( (status&0x0040)==0x0040) || !( (status&0x0004)==0x0004) ||  !((status&0x0080)==0x0080) ) return;
 203
 204     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 205     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 206     if (track->GetTPCNclsF()>0) {
 207       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 208     }
 209
 210     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 211
 212     Double_t mom=track->GetTPCmomentum();
 213
 214     for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 215       TH2 *h=(TH2*)fListQAtpc->At(ispecie);
 216       if (!h) continue;
 217       Double_t nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
 218       h->Fill(mom,nSigma);
 219     }
 220
 221     TH2 *h=(TH2*)fListQAtpc->At(AliPID::kSPECIES);
 222     if (h) {
 223       Double_t sig=track->GetTPCsignal();
 224       h->Fill(mom,sig);
 225     }
 226   }
 227 }
 228
 229 //______________________________________________________________________________
 230 void AliAnalysisTaskPIDqa::FillTOFqa()
 231 {
 232   AliVEvent *event=InputEvent();
 233
 234   Int_t ntracks=event->GetNumberOfTracks();
 235   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 236     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 237
 238     //
 239     //basic track cuts
 240     //
 241     ULong_t status=track->GetStatus();
 242     // not that nice. status bits not in virtual interface
 243     // TPC refit + ITS refit +
 244     // TOF out + TOFpid +
 245     // kTIME
 246     if (!((status&0x0040)==0x0040) || !((status&0x0004)==0x0004)
 247         || !((status&0x2000)==0x2000) || !((status&0x8000)==0x8000)
 248         || !((status&0x80000000)==0x80000000) ) return;
 249
 250     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 251     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 252     if (track->GetTPCNclsF()>0) {
 253       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 254     }
 255
 256     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 257
 258
 259     Double_t mom=track->P();
 260
 261     for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 262       TH2 *h=(TH2*)fListQAtof->At(ispecie);
 263       if (!h) continue;
 264       Double_t nSigma=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
 265       h->Fill(mom,nSigma);
 266     }
 267
 268     TH2 *h=(TH2*)fListQAtof->At(AliPID::kSPECIES);
 269     if (h) {
 270       Double_t sig=track->GetTOFsignal();
 271       h->Fill(mom,sig);
 272     }
 273
 274   }
 275 }
 276
 277 //______________________________________________________________________________
 278 void AliAnalysisTaskPIDqa::SetupITSqa()
 279 {
 280   //
 281   // Create the ITS qa objects
 282   //
 283
 284   TVectorD *vX=MakeLogBinning(200,.1,30);
 285
 286   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 287     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_ITS_%s",AliPID::ParticleName(ispecie)),
 288                               Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 289                               vX->GetNrows()-1,vX->GetMatrixArray(),
 290                               200,-10,10);
 291     fListQAits->Add(hNsigmaP);
 292   }
 293
 294
 295   TH2F *hSig = new TH2F("hSigP_ITS",
 296                         "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
 297                         vX->GetNrows()-1,vX->GetMatrixArray(),
 298                         300,0,300);
 299
 300   fListQAits->Add(hSig);
 301
 302 }
 303
 304 //______________________________________________________________________________
 305 void AliAnalysisTaskPIDqa::SetupTPCqa()
 306 {
 307   //
 308   // Create the TPC qa objects
 309   //
 310
 311   TVectorD *vX=MakeLogBinning(200,.1,30);
 312
 313   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 314     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_%s",AliPID::ParticleName(ispecie)),
 315                               Form("TPC n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 316                               vX->GetNrows()-1,vX->GetMatrixArray(),
 317                               200,-10,10);
 318     fListQAtpc->Add(hNsigmaP);
 319   }
 320
 321
 322   TH2F *hSig = new TH2F("hSigP_TPC",
 323                         "TPC signal vs. p;p [GeV]; TPC signal [arb. units]",
 324                         vX->GetNrows()-1,vX->GetMatrixArray(),
 325                         300,0,300);
 326
 327   fListQAtpc->Add(hSig);
 328
 329 }
 330
 331 //______________________________________________________________________________
 332 void AliAnalysisTaskPIDqa::SetupTRDqa()
 333 {
 334   //
 335   // Create the TRD qa objects
 336   //
 337
 338 }
 339
 340 //______________________________________________________________________________
 341 void AliAnalysisTaskPIDqa::SetupTOFqa()
 342 {
 343   //
 344   // Create the TOF qa objects
 345   //
 346
 347   TVectorD *vX=MakeLogBinning(200,.1,30);
 348
 349   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 350     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_%s",AliPID::ParticleName(ispecie)),
 351                               Form("TOF n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 352                               vX->GetNrows()-1,vX->GetMatrixArray(),
 353                               200,-10,10);
 354     fListQAtof->Add(hNsigmaP);
 355   }
 356
 357
 358   TH2F *hSig = new TH2F("hSigP_TOF",
 359                         "TOF signal vs. p;p [GeV]; TOF signal [arb. units]",
 360                         vX->GetNrows()-1,vX->GetMatrixArray(),
 361                         300,0,300);
 362
 363   fListQAtof->Add(hSig);
 364
 365 }
 366
 367 //______________________________________________________________________________
 368 TVectorD* AliAnalysisTaskPIDqa::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
 369 {
 370   //
 371   // Make logarithmic binning
 372   // the user has to delete the array afterwards!!!
 373   //
 374
 375   //check limits
 376   if (xmin<1e-20 || xmax<1e-20){
 377     AliError("For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
 378     return MakeLinBinning(nbinsX, xmin, xmax);
 379   }
 380   if (xmax<xmin){
 381     Double_t tmp=xmin;
 382     xmin=xmax;
 383     xmax=tmp;
 384   }
 385   TVectorD *binLim=new TVectorD(nbinsX+1);
 386   Double_t first=xmin;
 387   Double_t last=xmax;
 388   Double_t expMax=TMath::Log(last/first);
 389   for (Int_t i=0; i<nbinsX+1; ++i){
 390     (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
 391   }
 392   return binLim;
 393 }
 394
 395 //______________________________________________________________________________
 396 TVectorD* AliAnalysisTaskPIDqa::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
 397 {
 398   //
 399   // Make linear binning
 400   // the user has to delete the array afterwards!!!
 401   //
 402   if (xmax<xmin){
 403     Double_t tmp=xmin;
 404     xmin=xmax;
 405     xmax=tmp;
 406   }
 407   TVectorD *binLim=new TVectorD(nbinsX+1);
 408   Double_t first=xmin;
 409   Double_t last=xmax;
 410   Double_t binWidth=(last-first)/nbinsX;
 411   for (Int_t i=0; i<nbinsX+1; ++i){
 412     (*binLim)[i]=first+binWidth*(Double_t)i;
 413   }
 414   return binLim;
 415 }
 416
 417 //_____________________________________________________________________________
 418 TVectorD* AliAnalysisTaskPIDqa::MakeArbitraryBinning(const char* bins)
 419 {
 420   //
 421   // Make arbitrary binning, bins separated by a ','
 422   //
 423   TString limits(bins);
 424   if (limits.IsNull()){
 425     AliError("Bin Limit string is empty, cannot add the variable");
 426     return 0x0;
 427   }
 428
 429   TObjArray *arr=limits.Tokenize(",");
 430   Int_t nLimits=arr->GetEntries();
 431   if (nLimits<2){
 432     AliError("Need at leas 2 bin limits, cannot add the variable");
 433     delete arr;
 434     return 0x0;
 435   }
 436
 437   TVectorD *binLimits=new TVectorD(nLimits);
 438   for (Int_t iLim=0; iLim<nLimits; ++iLim){
 439     (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
 440   }
 441
 442   delete arr;
 443   return binLimits;
 444 }