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) AliError("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 & AliVTrack::kITSrefit)==AliVTrack::kITSrefit ) &&
 169            ( (status & AliVTrack::kITSpid  )==AliVTrack::kITSpid   ) )) continue;
 170     Double_t mom=track->P();
 171
 172     for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 173       TH2 *h=(TH2*)fListQAits->At(ispecie);
 174       if (!h) continue;
 175       Double_t nSigma=fPIDResponse->NumberOfSigmasITS(track, (AliPID::EParticleType)ispecie);
 176       h->Fill(mom,nSigma);
 177     }
 178
 179     TH2 *h=(TH2*)fListQAits->At(AliPID::kSPECIES);
 180     if (h) {
 181       Double_t sig=track->GetITSsignal();
 182       h->Fill(mom,sig);
 183     }
 184
 185   }
 186 }
 187
 188 //______________________________________________________________________________
 189 void AliAnalysisTaskPIDqa::FillTPCqa()
 190 {
 191   AliVEvent *event=InputEvent();
 192
 193   Int_t ntracks=event->GetNumberOfTracks();
 194   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 195     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 196
 197     //
 198     //basic track cuts
 199     //
 200     ULong_t status=track->GetStatus();
 201     // not that nice. status bits not in virtual interface
 202     // TPC refit + ITS refit + TPC pid
 203     if (!( (status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 204         !( (status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
 205         !( (status & AliVTrack::kTPCpid  ) == AliVTrack::kTPCpid  ) ) continue;
 206
 207     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 208     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 209     if (track->GetTPCNclsF()>0) {
 210       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 211     }
 212
 213     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 214
 215     Double_t mom=track->GetTPCmomentum();
 216
 217     for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 218       TH2 *h=(TH2*)fListQAtpc->At(ispecie);
 219       if (!h) continue;
 220       Double_t nSigma=fPIDResponse->NumberOfSigmasTPC(track, (AliPID::EParticleType)ispecie);
 221       h->Fill(mom,nSigma);
 222     }
 223
 224     TH2 *h=(TH2*)fListQAtpc->At(AliPID::kSPECIES);
 225     if (h) {
 226       Double_t sig=track->GetTPCsignal();
 227       h->Fill(mom,sig);
 228     }
 229   }
 230 }
 231
 232 //______________________________________________________________________________
 233 void AliAnalysisTaskPIDqa::FillTOFqa()
 234 {
 235   AliVEvent *event=InputEvent();
 236
 237   Int_t ntracks=event->GetNumberOfTracks();
 238   for(Int_t itrack = 0; itrack < ntracks; itrack++){
 239     AliVTrack *track=(AliVTrack*)event->GetTrack(itrack);
 240
 241     //
 242     //basic track cuts
 243     //
 244     ULong_t status=track->GetStatus();
 245     // not that nice. status bits not in virtual interface
 246     // TPC refit + ITS refit +
 247     // TOF out + TOFpid +
 248     // kTIME
 249     if (!((status & AliVTrack::kTPCrefit) == AliVTrack::kTPCrefit) ||
 250         !((status & AliVTrack::kITSrefit) == AliVTrack::kITSrefit) ||
 251         !((status & AliVTrack::kTOFout  ) == AliVTrack::kTOFout  ) ||
 252         !((status & AliVTrack::kTOFpid  ) == AliVTrack::kTOFpid  ) ||
 253         !((status & AliVTrack::kTIME    ) == AliVTrack::kTIME    ) ) continue;
 254
 255     Float_t nCrossedRowsTPC = track->GetTPCClusterInfo(2,1);
 256     Float_t  ratioCrossedRowsOverFindableClustersTPC = 1.0;
 257     if (track->GetTPCNclsF()>0) {
 258       ratioCrossedRowsOverFindableClustersTPC = nCrossedRowsTPC/track->GetTPCNclsF();
 259     }
 260
 261     if ( nCrossedRowsTPC<70 || ratioCrossedRowsOverFindableClustersTPC<.8 ) continue;
 262
 263
 264     Double_t mom=track->P();
 265
 266     for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 267       TH2 *h=(TH2*)fListQAtof->At(ispecie);
 268       if (!h) continue;
 269       Double_t nSigma=fPIDResponse->NumberOfSigmasTOF(track, (AliPID::EParticleType)ispecie);
 270       h->Fill(mom,nSigma);
 271     }
 272
 273     TH2 *h=(TH2*)fListQAtof->At(AliPID::kSPECIES);
 274     if (h) {
 275       Double_t sig=track->GetTOFsignal();
 276       h->Fill(mom,sig);
 277     }
 278
 279   }
 280 }
 281
 282 //______________________________________________________________________________
 283 void AliAnalysisTaskPIDqa::SetupITSqa()
 284 {
 285   //
 286   // Create the ITS qa objects
 287   //
 288
 289   TVectorD *vX=MakeLogBinning(200,.1,30);
 290
 291   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 292     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_ITS_%s",AliPID::ParticleName(ispecie)),
 293                               Form("ITS n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 294                               vX->GetNrows()-1,vX->GetMatrixArray(),
 295                               200,-10,10);
 296     fListQAits->Add(hNsigmaP);
 297   }
 298
 299
 300   TH2F *hSig = new TH2F("hSigP_ITS",
 301                         "ITS signal vs. p;p [GeV]; ITS signal [arb. units]",
 302                         vX->GetNrows()-1,vX->GetMatrixArray(),
 303                         300,0,300);
 304
 305   fListQAits->Add(hSig);
 306
 307 }
 308
 309 //______________________________________________________________________________
 310 void AliAnalysisTaskPIDqa::SetupTPCqa()
 311 {
 312   //
 313   // Create the TPC qa objects
 314   //
 315
 316   TVectorD *vX=MakeLogBinning(200,.1,30);
 317
 318   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 319     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TPC_%s",AliPID::ParticleName(ispecie)),
 320                               Form("TPC n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 321                               vX->GetNrows()-1,vX->GetMatrixArray(),
 322                               200,-10,10);
 323     fListQAtpc->Add(hNsigmaP);
 324   }
 325
 326
 327   TH2F *hSig = new TH2F("hSigP_TPC",
 328                         "TPC signal vs. p;p [GeV]; TPC signal [arb. units]",
 329                         vX->GetNrows()-1,vX->GetMatrixArray(),
 330                         300,0,300);
 331
 332   fListQAtpc->Add(hSig);
 333
 334 }
 335
 336 //______________________________________________________________________________
 337 void AliAnalysisTaskPIDqa::SetupTRDqa()
 338 {
 339   //
 340   // Create the TRD qa objects
 341   //
 342
 343 }
 344
 345 //______________________________________________________________________________
 346 void AliAnalysisTaskPIDqa::SetupTOFqa()
 347 {
 348   //
 349   // Create the TOF qa objects
 350   //
 351
 352   TVectorD *vX=MakeLogBinning(200,.1,30);
 353
 354   for (Int_t ispecie=0; ispecie<AliPID::kSPECIES; ++ispecie){
 355     TH2F *hNsigmaP = new TH2F(Form("hNsigmaP_TOF_%s",AliPID::ParticleName(ispecie)),
 356                               Form("TOF n#sigma %s vs. p;p [GeV]; n#sigma",AliPID::ParticleName(ispecie)),
 357                               vX->GetNrows()-1,vX->GetMatrixArray(),
 358                               200,-10,10);
 359     fListQAtof->Add(hNsigmaP);
 360   }
 361
 362
 363   TH2F *hSig = new TH2F("hSigP_TOF",
 364                         "TOF signal vs. p;p [GeV]; TOF signal [arb. units]",
 365                         vX->GetNrows()-1,vX->GetMatrixArray(),
 366                         300,0,300);
 367
 368   fListQAtof->Add(hSig);
 369
 370 }
 371
 372 //______________________________________________________________________________
 373 TVectorD* AliAnalysisTaskPIDqa::MakeLogBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
 374 {
 375   //
 376   // Make logarithmic binning
 377   // the user has to delete the array afterwards!!!
 378   //
 379
 380   //check limits
 381   if (xmin<1e-20 || xmax<1e-20){
 382     AliError("For Log binning xmin and xmax must be > 1e-20. Using linear binning instead!");
 383     return MakeLinBinning(nbinsX, xmin, xmax);
 384   }
 385   if (xmax<xmin){
 386     Double_t tmp=xmin;
 387     xmin=xmax;
 388     xmax=tmp;
 389   }
 390   TVectorD *binLim=new TVectorD(nbinsX+1);
 391   Double_t first=xmin;
 392   Double_t last=xmax;
 393   Double_t expMax=TMath::Log(last/first);
 394   for (Int_t i=0; i<nbinsX+1; ++i){
 395     (*binLim)[i]=first*TMath::Exp(expMax/nbinsX*(Double_t)i);
 396   }
 397   return binLim;
 398 }
 399
 400 //______________________________________________________________________________
 401 TVectorD* AliAnalysisTaskPIDqa::MakeLinBinning(Int_t nbinsX, Double_t xmin, Double_t xmax)
 402 {
 403   //
 404   // Make linear binning
 405   // the user has to delete the array afterwards!!!
 406   //
 407   if (xmax<xmin){
 408     Double_t tmp=xmin;
 409     xmin=xmax;
 410     xmax=tmp;
 411   }
 412   TVectorD *binLim=new TVectorD(nbinsX+1);
 413   Double_t first=xmin;
 414   Double_t last=xmax;
 415   Double_t binWidth=(last-first)/nbinsX;
 416   for (Int_t i=0; i<nbinsX+1; ++i){
 417     (*binLim)[i]=first+binWidth*(Double_t)i;
 418   }
 419   return binLim;
 420 }
 421
 422 //_____________________________________________________________________________
 423 TVectorD* AliAnalysisTaskPIDqa::MakeArbitraryBinning(const char* bins)
 424 {
 425   //
 426   // Make arbitrary binning, bins separated by a ','
 427   //
 428   TString limits(bins);
 429   if (limits.IsNull()){
 430     AliError("Bin Limit string is empty, cannot add the variable");
 431     return 0x0;
 432   }
 433
 434   TObjArray *arr=limits.Tokenize(",");
 435   Int_t nLimits=arr->GetEntries();
 436   if (nLimits<2){
 437     AliError("Need at leas 2 bin limits, cannot add the variable");
 438     delete arr;
 439     return 0x0;
 440   }
 441
 442   TVectorD *binLimits=new TVectorD(nLimits);
 443   for (Int_t iLim=0; iLim<nLimits; ++iLim){
 444     (*binLimits)[iLim]=(static_cast<TObjString*>(arr->At(iLim)))->GetString().Atof();
 445   }
 446
 447   delete arr;
 448   return binLimits;
 449 }