PWG1/AliTPCtaskPID.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 //
  17 //
  18 //Task for analysis if TPC dEdx and PID information
  19 //
  20 //-----------------------------------------------------------------------
  21 // Author : M.Ivanov  marian.ivanov@cern.ch -
  22 //-----------------------------------------------------------------------
  23
  24
  25 // 3 6D histograms  - THnSparse created in the task:
  26 // TPC raw dEdx
  27 // TPC normalized dEdx (dEdx_rec/dNdx_mc)
  28 // TPC PID probabilities
  29 //
  30 // The values are binned in following variables:
  31 // Some of them are correlated - but THnSpase handle it
  32 //                               ~ 14 MBy per object needed
  33 //
  34 // 0 - MC particle species as defined in the AliPID - negatives+5 <0,9>
  35 // 1 - momenta - at the entrance of the TPC
  36 // 2 - tan lambda- fP[3]
  37 // 3 - betagamma
  38 // 4 - npoints
  39 // 5 - measurement - dEdx, dEdx/BB resp.  PID probability
  40 // 6 - BB resp. rec particle
  41
  42
  43 //
  44 // ROOT includes
  45 #include <TChain.h>
  46 #include <TMath.h>
  47 #include <TVectorD.h>
  48 #include <TSystem.h>
  49 #include <TFile.h>
  50 #include <TParticle.h>
  51
  52 // ALIROOT includes
  53 #include <AliAnalysisManager.h>
  54 #include <AliESDInputHandler.h>
  55 #include "AliStack.h"
  56 #include "AliMCEvent.h"
  57 #include "AliMCEventHandler.h"
  58 #include "AliMathBase.h"
  59
  60 #include <AliESD.h>
  61 #include "AliExternalTrackParam.h"
  62 #include "AliTracker.h"
  63 #include "AliTPCseed.h"
  64 //
  65 #include "AliTPCtaskPID.h"
  66 //
  67 #include <THnSparse.h>
  68
  69 //
  70
  71 // STL includes
  72 #include <iostream>
  73
  74 using namespace std;
  75
  76
  77 ClassImp(AliTPCtaskPID)
  78
  79 //________________________________________________________________________
  80 AliTPCtaskPID::AliTPCtaskPID() :
  81   AliAnalysisTask(),
  82   fMCinfo(0),     //! MC event handler
  83   fESD(0),
  84   fList(0),
  85   fTPCsignal(0),
  86   fTPCsignalNorm(0),
  87   fTPCr(0)
  88 {
  89   //
  90   // Default constructor (should not be used)
  91   //
  92 }
  93
  94 AliTPCtaskPID::AliTPCtaskPID(const AliTPCtaskPID& info) :
  95   AliAnalysisTask(info),
  96   fMCinfo(info.fMCinfo),     //! MC event handler
  97   fESD(info.fESD),        //!
  98   fList(0),
  99   fTPCsignal(0),
 100   fTPCsignalNorm(0),
 101   fTPCr(0)
 102 {
 103   //
 104   // Dummy Copy  constructor - no copy constructor for THnSparse
 105   //
 106   fList = (TObjArray*)(info.fList->Clone());
 107 }
 108
 109
 110
 111 //________________________________________________________________________
 112 AliTPCtaskPID::AliTPCtaskPID(const char *name) :
 113   AliAnalysisTask(name, "AliTPCtaskPID"),
 114   fMCinfo(0),     //! MC event handler
 115   fESD(0),
 116   fList(0),
 117   fTPCsignal(0),
 118   fTPCsignalNorm(0),
 119   fTPCr(0)
 120 {
 121   //
 122   // Normal constructor
 123   //
 124   // Input slot #0 works with a TChain
 125   DefineInput(0, TChain::Class());
 126   // Output slot #0 writes into a TList
 127   DefineOutput(0, TObjArray::Class());
 128   //
 129   //make histos
 130   Init();
 131 }
 132
 133 void AliTPCtaskPID::Init(){
 134   //
 135   // Init dEdx histogram
 136   // Dimensions
 137   // 0 - particle specie as defined in the AliPID - negatives+5 <0,9>
 138   // 1 - momenta - at the entrance of the TPC
 139   // 2 - tan lambda- fP[3]
 140   // 3 - betagamma
 141   // 4 - npoints
 142   // 5 - measurement - dEdx or dEdx/BB
 143   // 6 - BB
 144   //
 145   Double_t xmin[7],  xmax[7];
 146   Int_t    nbins[7];
 147   // pid
 148   nbins[0]=10;
 149   xmin[0]=0; xmax[0]=10;
 150   // momenta
 151   nbins[1]=50;
 152   xmin[1]=0.1; xmax[1]=100;
 153   //pseudorapidity
 154   nbins[2]=40;
 155   xmin[2]=-1.4; xmax[2]=1.4;
 156   //
 157   // betagamma
 158   //
 159   nbins[3]=100;
 160   xmin[3]=0.1; xmax[3]=1000;
 161   //
 162   nbins[4]=11;
 163   xmin[4] =50; xmax[4]=160;
 164   //
 165   //
 166   nbins[6]=60;
 167   xmin[6]=1; xmax[6]=6;
 168
 169   nbins[5]=400;
 170   xmin[5]=20; xmax[5]=400;
 171   fTPCsignal = new THnSparseF("TPC signal","TPC signal",7,nbins,xmin,xmax);
 172   nbins[5]=100;
 173   xmin[5]=25; xmax[5]=75;
 174   fTPCsignalNorm = new THnSparseF("TPC signal Norm","TPC signal Norm",7,nbins,xmin,xmax);
 175   //
 176   nbins[5]=256;
 177   xmin[5]=-0.1; xmax[5]=1.1;
 178   nbins[6]=10;
 179   xmin[6]=0; xmax[6]=10;
 180   fTPCr = new THnSparseF("TPC pid probability ","TPC pid probability",7,nbins,xmin,xmax);
 181   //
 182   //
 183   BinLogX(fTPCsignal->GetAxis(1));
 184   BinLogX(fTPCsignal->GetAxis(3));
 185   BinLogX(fTPCsignal->GetAxis(6));
 186   BinLogX(fTPCsignalNorm->GetAxis(1));
 187   BinLogX(fTPCsignalNorm->GetAxis(3));
 188   BinLogX(fTPCsignalNorm->GetAxis(6));
 189   BinLogX(fTPCr->GetAxis(1));
 190   BinLogX(fTPCr->GetAxis(3));
 191
 192   char *hisAxisName[7] ={"pid","p (GeV/c)","#eta","#beta#gamma","Number of cluster","dEdx_{rec},dedx_{mc}"};
 193   char *hisAxisNameNorm[7] ={"pid","p (GeV/c)","#eta","#beta#gamma","Number of cluster","dEdx_{rec}/dEdx_{mc},dedx_{mc}"};
 194   char *hisAxisNameR[7] ={"pid","p (GeV/c)","#eta","#beta#gamma","Number of cluster","TPCr","pid2"};
 195   //
 196   for (Int_t i=0;i<7;i++) {
 197     fTPCsignal->GetAxis(i)->SetTitle(hisAxisName[i]);
 198     fTPCsignal->GetAxis(i)->SetName(hisAxisName[i]);
 199     fTPCsignalNorm->GetAxis(i)->SetTitle(hisAxisNameNorm[i]);
 200     fTPCsignalNorm->GetAxis(i)->SetName(hisAxisNameNorm[i]);
 201     fTPCr->GetAxis(i)->SetTitle(hisAxisNameR[i]);
 202     fTPCr->GetAxis(i)->SetName(hisAxisNameR[i]);
 203   }
 204
 205   fList = new TObjArray(3);
 206   fList->AddAt(fTPCsignal,0);
 207   fList->AddAt(fTPCsignalNorm,1);
 208   fList->AddAt(fTPCr,2);
 209 }
 210
 211
 212
 213
 214 AliTPCtaskPID::~AliTPCtaskPID(){
 215   //
 216   //
 217   //
 218   delete fTPCsignal;
 219   delete fTPCsignalNorm;
 220   delete fTPCr;
 221 }
 222
 223
 224 //________________________________________________________________________
 225 void AliTPCtaskPID::ConnectInputData(Option_t *)
 226 {
 227   //
 228   // Connect the input data
 229   //
 230
 231   TTree* tree=dynamic_cast<TTree*>(GetInputData(0));
 232   if (!tree) {
 233     //Printf("ERROR: Could not read chain from input slot 0");
 234   }
 235   else {
 236     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 237     if (!esdH) {
 238       //Printf("ERROR: Could not get ESDInputHandler");
 239     }
 240     else {
 241       fESD = esdH->GetEvent();
 242       //Printf("*** CONNECTED NEW EVENT ****");
 243     }
 244   }
 245   AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 246   mcinfo->SetReadTR(kTRUE);
 247
 248   fMCinfo = mcinfo->MCEvent();
 249
 250
 251 }
 252
 253
 254
 255
 256
 257
 258 //________________________________________________________________________
 259 void AliTPCtaskPID::CreateOutputObjects()
 260 {
 261   //
 262   // Connect the output objects
 263   //
 264
 265 }
 266
 267
 268 //________________________________________________________________________
 269 void AliTPCtaskPID::Exec(Option_t *) {
 270   //
 271   // Execute analysis for current event
 272   //
 273
 274
 275   // If MC has been connected
 276
 277   if (!fMCinfo){
 278     cout << "Not MC info\n" << endl;
 279   }else{
 280     ProcessMCInfo();
 281     //mcinfo->Print();
 282     //DumpInfo();
 283   }
 284   //
 285   PostData(0, fList);
 286 }
 287
 288
 289
 290
 291 //________________________________________________________________________
 292 void AliTPCtaskPID::Terminate(Option_t *) {
 293     //
 294     // Terminate loop
 295     //
 296   //
 297   return;
 298 }
 299
 300
 301
 302
 303
 304
 305
 306 void  AliTPCtaskPID::ProcessMCInfo(){
 307   //
 308   //
 309   //
 310   //
 311
 312   if (!fTPCsignal) Init();
 313   Int_t npart   = fMCinfo->GetNumberOfTracks();
 314   Int_t ntracks = fESD->GetNumberOfTracks();
 315   if (npart<=0) return;
 316   if (ntracks<=0) return;
 317   //
 318   //
 319   TParticle * particle= new TParticle;
 320   TClonesArray * trefs = new TClonesArray("AliTrackReference");
 321
 322   for (Int_t itrack=0;itrack<ntracks;itrack++){
 323     AliESDtrack *track = fESD->GetTrack(itrack);
 324     const AliExternalTrackParam *in=track->GetInnerParam();
 325     const AliExternalTrackParam *out=track->GetOuterParam();
 326     if (!in) continue;
 327     Int_t ipart = TMath::Abs(track->GetLabel());
 328     //
 329     Int_t status = fMCinfo->GetParticleAndTR(ipart, particle, trefs);
 330     if (status<0) continue;
 331     if (!particle) continue;
 332     if (!trefs) continue;
 333     //
 334     //
 335     Double_t mom = in->GetP();
 336     if (track->GetP()>5)  mom= track->GetP();
 337     if (out&&out->GetX()<300)  mom= (in->GetP()+out->GetP())*0.5;
 338     Double_t dedx=track->GetTPCsignal();
 339     Double_t mass = particle->GetMass();
 340     Double_t bg  =mom/mass;
 341     Double_t betheBloch = AliMathBase::BetheBlochAleph(bg);
 342     //
 343     // Fill histos
 344     //
 345     Double_t x[5];
 346     //PID
 347     Int_t pdg = particle->GetPdgCode();
 348     for (Int_t iType=0;iType<5;iType++) {
 349       if (AliPID::ParticleCode(iType)==TMath::Abs(pdg)){
 350         x[0]=iType;
 351         if (TDatabasePDG::Instance()->GetParticle(pdg)->Charge()<0) x[0]+=5;
 352       }
 353     }
 354     x[1]= mom;
 355     x[2]= -0.5*TMath::Log((track->P()+track->Pz())/(track->P()-track->Pz()));
 356     x[3]= bg;
 357     x[4]=track->GetTPCsignalN();
 358     x[5]= dedx;
 359     x[6]= betheBloch;
 360     //
 361     fTPCsignal->Fill(x);
 362     //
 363     x[5]=dedx/betheBloch;
 364     fTPCsignalNorm->Fill(x);
 365     for (Int_t ipart2=0;ipart2<5;ipart2++){
 366       x[6]=ipart2;
 367       Double_t tpcpid[AliPID::kSPECIES];
 368       track->GetTPCpid(tpcpid);
 369       x[5]=tpcpid[ipart2];
 370       fTPCr->Fill(x);
 371     }
 372   }
 373 }
 374
 375
 376
 377
 378 void AliTPCtaskPID::FinishTaskOutput()
 379 {
 380   //
 381   // According description in AliAnalisysTask this method is call
 382   // on the slaves before sending data
 383   //
 384   Terminate("slave");
 385   gSystem->Exec("pwd");
 386 }
 387
 388
 389
 390
 391 void AliTPCtaskPID::BinLogX(TAxis *axis) {
 392   //
 393   //
 394   //
 395   Int_t bins = axis->GetNbins();
 396
 397   Double_t from = axis->GetXmin();
 398   Double_t to   = axis->GetXmax();
 399   Double_t *new_bins = new Double_t[bins + 1];
 400
 401   new_bins[0] = from;
 402   Double_t factor = pow(to/from, 1./bins);
 403
 404   for (Int_t i = 1; i <= bins; i++) {
 405     new_bins[i] = factor * new_bins[i-1];
 406   }
 407   axis->Set(bins, new_bins);
 408   delete new_bins;
 409 }