PWGPP/TPC/AliTPCComparisonPID.cxx

   1 //
   2 //
   3 //
   4
   5 //
   6 // ROOT includes
   7 #include <TChain.h>
   8 #include <TMath.h>
   9 #include <TVectorD.h>
  10 #include <TSystem.h>
  11 #include <TFile.h>
  12 #include <TParticle.h>
  13
  14 // ALIROOT includes
  15 #include <TTreeStream.h>
  16 #include <AliAnalysisManager.h>
  17 #include <AliESDInputHandler.h>
  18 #include "AliStack.h"
  19 #include "AliMCEvent.h"
  20 #include "AliMCEventHandler.h"
  21 #include "AliMathBase.h"
  22
  23 #include <AliESD.h>
  24 #include "AliExternalTrackParam.h"
  25 #include "AliTracker.h"
  26 #include "AliTPCseed.h"
  27 //
  28 #include "AliTPCComparisonPID.h"
  29 //
  30 #include <THnSparse.h>
  31
  32 //
  33
  34 // STL includes
  35 #include <iostream>
  36
  37 using namespace std;
  38
  39 ClassImp(AliTPCComparisonPID)
  40
  41 //________________________________________________________________________
  42 AliTPCComparisonPID::AliTPCComparisonPID() :
  43   AliAnalysisTask(),
  44   fMCinfo(0),     //! MC event handler
  45   fESD(0),
  46   fTPCsignal(0),
  47   fTPCsignalNorm(0),
  48   fDebugStreamer(0),
  49   fStreamLevel(0),
  50   fDebugLevel(0),
  51   fDebugOutputPath()
  52 {
  53   //
  54   // Default constructor (should not be used)
  55   //
  56 }
  57
  58 AliTPCComparisonPID::AliTPCComparisonPID(const AliTPCComparisonPID& info) :
  59   AliAnalysisTask(info),
  60   fMCinfo(info.fMCinfo),     //! MC event handler
  61   fESD(info.fESD),        //!
  62   fTPCsignal(0),
  63   fTPCsignalNorm(0),
  64   //
  65   fDebugStreamer(0),
  66   fStreamLevel(0),
  67   fDebugLevel(),
  68   fDebugOutputPath()
  69 {
  70   //
  71   // Dummy Copy  constructor - no copy constructor for THnSparse
  72   //
  73 }
  74
  75
  76
  77 //________________________________________________________________________
  78 AliTPCComparisonPID::AliTPCComparisonPID(const char *name) :
  79   AliAnalysisTask(name, "AliTPCComparisonPID"),
  80   fMCinfo(0),     //! MC event handler
  81   fESD(0),
  82   fTPCsignal(0),
  83   fTPCsignalNorm(0),
  84   fDebugStreamer(0),
  85   fStreamLevel(0),
  86   fDebugLevel(0),
  87   fDebugOutputPath()
  88 {
  89   //
  90   // Normal constructor
  91   //
  92   // Input slot #0 works with a TChain
  93   DefineInput(0, TChain::Class());
  94   // Output slot #0 writes into a TList
  95   DefineOutput(0, AliTPCComparisonPID::Class());
  96   //
  97   //make histos
  98   Init();
  99 }
 100
 101 void AliTPCComparisonPID::Init(){
 102   //
 103   // Init dEdx histogram
 104   // Dimensions
 105   // 0 - particle specie as defined in the AliPID - negatives+5 <0,9>
 106   // 1 - momenta - at the entrance of the TPC
 107   // 2 - tan lambda- fP[3]
 108   // 3 - betagamma
 109   // 4 - measurement - dEdx or dEdx/BB
 110   //
 111   Double_t xmin[5],  xmax[5];
 112   Int_t    nbins[5];
 113   // pid
 114   nbins[0]=10;
 115   xmin[0]=0; xmax[0]=10;
 116   // momenta
 117   nbins[1]=30;
 118   xmin[1]=0.1; xmax[1]=3;
 119   //P3
 120   nbins[2]=20;
 121   xmin[2]=-1.5; xmax[2]=1.5;
 122   //
 123   // log (betagamma)
 124   //
 125   nbins[3]=50;
 126   xmin[3]=0.1; xmax[3]=100;
 127   //
 128   //
 129   nbins[4]=400;
 130   xmin[4]=20; xmax[4]=400;
 131   fTPCsignal = new THnSparseF("TPC signal","TPC signal",5,nbins,xmin,xmax);
 132   nbins[4]=100;
 133   xmin[4]=25; xmax[4]=75;
 134   fTPCsignal = new THnSparseF("TPC signal Norm","TPC signal Norm",5,nbins,xmin,xmax);
 135   //
 136 }
 137
 138
 139
 140
 141
 142 AliTPCComparisonPID::~AliTPCComparisonPID(){
 143   //
 144   //
 145   //
 146   if (fDebugLevel>0)  printf("AliTPCComparisonPID::~AliTPCComparisonPID\n");
 147   if (fDebugStreamer) delete fDebugStreamer;
 148   fDebugStreamer=0;
 149   delete fTPCsignal;
 150   delete fTPCsignalNorm;
 151 }
 152
 153
 154 //________________________________________________________________________
 155 void AliTPCComparisonPID::ConnectInputData(Option_t *)
 156 {
 157   //
 158   // Connect the input data
 159   //
 160   if(fDebugLevel>3)
 161     cout << "AnalysisTaskTPCCluster::ConnectInputData()" << endl;
 162
 163   TTree* tree=dynamic_cast<TTree*>(GetInputData(0));
 164   if (!tree) {
 165     //Printf("ERROR: Could not read chain from input slot 0");
 166   }
 167   else {
 168     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 169     if (!esdH) {
 170       //Printf("ERROR: Could not get ESDInputHandler");
 171     }
 172     else {
 173       fESD = esdH->GetEvent();
 174       //Printf("*** CONNECTED NEW EVENT ****");
 175     }
 176   }
 177   AliMCEventHandler* mcinfo = (AliMCEventHandler*) (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
 178   mcinfo->SetReadTR(kTRUE);
 179
 180   fMCinfo = mcinfo->MCEvent();
 181
 182
 183 }
 184
 185
 186
 187
 188
 189
 190 //________________________________________________________________________
 191 void AliTPCComparisonPID::CreateOutputObjects()
 192 {
 193   //
 194   // Connect the output objects
 195   //
 196   if(fDebugLevel>3)
 197     cout << "AnalysisTaskTPCCluster::CreateOutputObjects()" << endl;
 198
 199 }
 200
 201
 202 //________________________________________________________________________
 203 void AliTPCComparisonPID::Exec(Option_t *) {
 204   //
 205   // Execute analysis for current event
 206   //
 207
 208   if(fDebugLevel>3)
 209     cout << "AliTPCComparisonPID::Exec()" << endl;
 210
 211
 212   // If MC has been connected
 213
 214   if (!fMCinfo){
 215     cout << "Not MC info\n" << endl;
 216   }else{
 217     ProcessMCInfo();
 218     //mcinfo->Print();
 219     //DumpInfo();
 220   }
 221   //
 222   PostData(0, this);
 223 }
 224
 225
 226
 227
 228 //________________________________________________________________________
 229 void AliTPCComparisonPID::Terminate(Option_t *) {
 230     //
 231     // Terminate loop
 232     //
 233   if(fDebugLevel>3)
 234     printf("AliTPCComparisonPID: Terminate() \n");
 235   //
 236   if (fDebugLevel>0) printf("AliMCtrackingTestTask::Terminate\n");
 237   if (fDebugStreamer) delete fDebugStreamer;
 238   fDebugStreamer = 0;
 239   return;
 240 }
 241
 242
 243
 244 TTreeSRedirector *AliTPCComparisonPID::GetDebugStreamer(){
 245   //
 246   // Get Debug streamer
 247   // In case debug streamer not yet initialized and StreamLevel>0 create new one
 248   //
 249   if (fStreamLevel==0) return 0;
 250   if (fDebugStreamer) return fDebugStreamer;
 251   TString dsName;
 252   dsName=GetName();
 253   dsName+="Debug.root";
 254   dsName.ReplaceAll(" ","");
 255   fDebugStreamer = new TTreeSRedirector(dsName.Data());
 256   return fDebugStreamer;
 257 }
 258
 259
 260
 261
 262 void  AliTPCComparisonPID::ProcessMCInfo(){
 263   //
 264   //
 265   //
 266   //
 267   Int_t npart   = fMCinfo->GetNumberOfTracks();
 268   Int_t ntracks = fESD->GetNumberOfTracks();
 269   if (npart<=0) return;
 270   if (ntracks<=0) return;
 271   //
 272   //
 273   TParticle * particle= new TParticle;
 274   TClonesArray * trefs = new TClonesArray("AliTrackReference");
 275
 276   for (Int_t itrack=0;itrack<ntracks;itrack++){
 277     AliESDtrack *track = fESD->GetTrack(itrack);
 278     const AliExternalTrackParam *in=track->GetInnerParam();
 279     if (!in) continue;
 280     Int_t ipart = TMath::Abs(track->GetLabel());
 281     //
 282     Int_t status = fMCinfo->GetParticleAndTR(ipart, particle, trefs);
 283     if (status<0) continue;
 284     if (!particle) continue;
 285     if (!trefs) continue;
 286     //
 287     //
 288     Double_t mom = in->GetP();
 289     Double_t dedx=track->GetTPCsignal();
 290     Double_t mass = particle->GetMass();
 291     Double_t bg  =mom/mass;
 292     Double_t betheBloch = AliMathBase::BetheBlochAleph(bg);
 293     //
 294     // Fill histos
 295     //
 296     Double_t x[5];
 297     //PID
 298     Int_t pdg = particle->GetPdgCode();
 299     for (Int_t iType=0;iType<5;iType++) {
 300       if (AliPID::ParticleCode(iType)==TMath::Abs(pdg)){
 301         x[0]=iType;
 302         if (pdg<0) x[0]+=5;
 303       }
 304     }
 305     x[1]= mom;
 306     x[2]= track->GetTgl();
 307     x[3]= TMath::Log(bg);
 308     x[4]= dedx;
 309     fTPCsignal->Fill(x);
 310     x[4]=dedx/betheBloch;
 311     fTPCsignalNorm->Fill(x);
 312   }
 313 }
 314
 315
 316
 317
 318 void AliTPCComparisonPID::FinishTaskOutput()
 319 {
 320   //
 321   // According description in AliAnalisysTask this method is call
 322   // on the slaves before sending data
 323   //
 324   Terminate("slave");
 325   gSystem->Exec("pwd");
 326   RegisterDebugOutput();
 327
 328 }
 329
 330
 331 void AliTPCComparisonPID::RegisterDebugOutput(){
 332   //
 333   //
 334   //
 335   //
 336   // store  - copy debug output to the destination position
 337   // currently ONLY for local copy
 338   TString dsName;
 339   dsName=GetName();
 340   dsName+="Debug.root";
 341   dsName.ReplaceAll(" ","");
 342   TString dsName2=fDebugOutputPath.Data();
 343   gSystem->MakeDirectory(dsName2.Data());
 344   dsName2+=gSystem->HostName();
 345   gSystem->MakeDirectory(dsName2.Data());
 346   dsName2+="/";
 347   dsName2+=gSystem->BaseName(gSystem->pwd());
 348   dsName2+="/";
 349   gSystem->MakeDirectory(dsName2.Data());
 350   dsName2+=dsName;
 351   AliInfo(Form("copy %s\t%s\n",dsName.Data(),dsName2.Data()));
 352   printf("copy %s\t%s\n",dsName.Data(),dsName2.Data());
 353   TFile::Cp(dsName.Data(),dsName2.Data());
 354 }