PWG1/AliPerformanceTPC.cxx

   1 //------------------------------------------------------------------------------
   2 // Implementation of AliPerformanceTPC class. It keeps information from
   3 // comparison of reconstructed and MC particle tracks. In addtion,
   4 // it keeps selection cuts used during comparison. The comparison
   5 // information is stored in the ROOT histograms. Analysis of these
   6 // histograms can be done by using Analyse() class function. The result of
   7 // the analysis (histograms/graphs) are stored in the folder which is
   8 // a data member of AliPerformanceTPC.
   9 //
  10 // Author: J.Otwinowski 04/02/2008
  11 //------------------------------------------------------------------------------
  12
  13 /*
  14
  15   // after running comparison task, read the file, and get component
  16   gROOT->LoadMacro("$ALICE_ROOT/PWG1/Macros/LoadMyLibs.C");
  17   LoadMyLibs();
  18
  19   TFile f("Output.root");
  20   AliPerformanceTPC * compObj = (AliPerformanceTPC*)coutput->FindObject("AliPerformanceTPC");
  21
  22   // analyse comparison data
  23   compObj->Analyse();
  24
  25   // the output histograms/graphs will be stored in the folder "folderTPC"
  26   compObj->GetAnalysisFolder()->ls("*");
  27
  28   // user can save whole comparison object (or only folder with anlysed histograms)
  29   // in the seperate output file (e.g.)
  30   TFile fout("Analysed_TPC.root","recreate");
  31   compObj->Write(); // compObj->GetAnalysisFolder()->Write();
  32   fout.Close();
  33
  34 */
  35
  36 #include "TCanvas.h"
  37 #include "TH1.h"
  38 #include "TH2.h"
  39 #include "TAxis.h"
  40 #include "TPostScript.h"
  41
  42 #include "AliPerformanceTPC.h"
  43 #include "AliESDEvent.h"
  44 #include "AliESDVertex.h"
  45 #include "AliESDtrack.h"
  46 #include "AliLog.h"
  47 #include "AliMCEvent.h"
  48 #include "AliHeader.h"
  49 #include "AliGenEventHeader.h"
  50 #include "AliStack.h"
  51 #include "AliMCInfoCuts.h"
  52 #include "AliRecInfoCuts.h"
  53 #include "AliTracker.h"
  54 #include "AliTreeDraw.h"
  55
  56 using namespace std;
  57
  58 ClassImp(AliPerformanceTPC)
  59
  60 //_____________________________________________________________________________
  61 AliPerformanceTPC::AliPerformanceTPC():
  62   AliPerformanceObject("AliPerformanceTPC"),
  63   fTPCHisto(0),
  64
  65   // Cuts
  66   fCutsRC(0),
  67   fCutsMC(0),
  68
  69   // histogram folder
  70   fAnalysisFolder(0)
  71 {
  72   Init();
  73 }
  74
  75 //_____________________________________________________________________________
  76 AliPerformanceTPC::AliPerformanceTPC(Char_t* name="AliPerformanceTPC", Char_t* title="AliPerformanceTPC",Int_t analysisMode=0,Bool_t hptGenerator=kFALSE):
  77   AliPerformanceObject(name,title),
  78   fTPCHisto(0),
  79
  80   // Cuts
  81   fCutsRC(0),
  82   fCutsMC(0),
  83
  84   // histogram folder
  85   fAnalysisFolder(0)
  86 {
  87   // named constructor
  88   //
  89   SetAnalysisMode(analysisMode);
  90   SetHptGenerator(hptGenerator);
  91
  92   Init();
  93 }
  94
  95 //_____________________________________________________________________________
  96 AliPerformanceTPC::~AliPerformanceTPC()
  97 {
  98   // destructor
  99
 100   if(fTPCHisto) delete fTPCHisto; fTPCHisto=0;
 101   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
 102 }
 103
 104 //_____________________________________________________________________________
 105 void AliPerformanceTPC::Init(){
 106   //
 107   // histogram bining
 108   //
 109
 110   // set pt bins
 111   Int_t nPtBins = 50;
 112   Double_t ptMin = 1.e-2, ptMax = 10.;
 113
 114   Double_t *binsPt = 0;
 115   if (IsHptGenerator())  {
 116     nPtBins = 100; ptMax = 100.;
 117     binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
 118   } else {
 119     binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
 120   }
 121
 122   /*
 123   Int_t nPtBins = 31;
 124   Double_t binsPt[32] = {0.,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8,2.0,2.25,2.5,2.75,3.,3.5,4.,5.,6.,8.,10.};
 125   Double_t ptMin = 0., ptMax = 10.;
 126
 127   if(IsHptGenerator() == kTRUE) {
 128     nPtBins = 100;
 129     ptMin = 0.; ptMax = 100.;
 130   }
 131   */
 132
 133   // nTPCClust:chi2PerTPCClust:nTPCClustFindRatio:eta:phi:pt
 134   Int_t binsTPCHisto[6]={160,100,100,30,144,nPtBins};
 135   Double_t minTPCHisto[6]={0., 0., 0., -1.5, 0., ptMin};
 136   Double_t maxTPCHisto[6]={160.,10.,1.2, 1.5, 2.*TMath::Pi(), ptMax};
 137
 138   fTPCHisto = new THnSparseF("fTPCHisto","nClust:chi2PerClust:nClust/nFindableClust:eta:phi:pt",6,binsTPCHisto,minTPCHisto,maxTPCHisto);
 139   fTPCHisto->SetBinEdges(5,binsPt);
 140
 141   fTPCHisto->GetAxis(0)->SetTitle("nClust");
 142   fTPCHisto->GetAxis(1)->SetTitle("chi2PerClust");
 143   fTPCHisto->GetAxis(2)->SetTitle("nClust/nFindableClust");
 144   fTPCHisto->GetAxis(3)->SetTitle("#eta");
 145   fTPCHisto->GetAxis(4)->SetTitle("#phi (rad)");
 146   fTPCHisto->GetAxis(5)->SetTitle("p_{T} (GeV/c)");
 147   fTPCHisto->Sumw2();
 148
 149   // Init cuts
 150   if(!fCutsMC)
 151     AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
 152   if(!fCutsRC)
 153     AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
 154
 155   // init folder
 156   fAnalysisFolder = CreateFolder("folderTPC","Analysis Resolution Folder");
 157 }
 158
 159 //_____________________________________________________________________________
 160 void AliPerformanceTPC::ProcessTPC(AliStack* const stack, AliESDtrack *const esdTrack)
 161 {
 162   if(!esdTrack) return;
 163
 164   // Fill TPC only resolution comparison information
 165   const AliExternalTrackParam *track = esdTrack->GetTPCInnerParam();
 166   if(!track) return;
 167
 168   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 169   esdTrack->GetImpactParametersTPC(dca,cov);
 170
 171   //Float_t q = esdTrack->Charge();
 172   Float_t pt = track->Pt();
 173   Float_t eta = track->Eta();
 174   Float_t phi = track->Phi();
 175   Int_t nClust = esdTrack->GetTPCclusters(0);
 176   Int_t nFindableClust = esdTrack->GetTPCNclsF();
 177
 178   Float_t chi2PerCluster = 0.;
 179   if(nClust>0.) chi2PerCluster = esdTrack->GetTPCchi2()/Float_t(nClust);
 180
 181   Float_t clustPerFindClust = 0.;
 182   if(nFindableClust>0.) clustPerFindClust = Float_t(nClust)/nFindableClust;
 183
 184   if(TMath::Abs(dca[0])<fCutsRC->GetMaxDCAToVertexXY() && TMath::Abs(dca[1])<fCutsRC->GetMaxDCAToVertexZ()) {
 185     Double_t vTPCHisto[6] = {nClust,chi2PerCluster,clustPerFindClust,eta,phi,pt};
 186     fTPCHisto->Fill(vTPCHisto);
 187   }
 188
 189   //
 190   // Fill rec vs MC information
 191   //
 192   if(!stack) return;
 193
 194 }
 195
 196 //_____________________________________________________________________________
 197 void AliPerformanceTPC::ProcessTPCITS(AliStack* const /*stack*/, AliESDtrack *const /*esdTrack*/)
 198 {
 199   // Fill comparison information (TPC+ITS)
 200   AliDebug(AliLog::kWarning, "Warning: Not implemented");
 201 }
 202
 203 //_____________________________________________________________________________
 204 void AliPerformanceTPC::ProcessConstrained(AliStack* const /*stack*/, AliESDtrack *const /*esdTrack*/)
 205 {
 206   // Fill comparison information (constarained parameters)
 207   AliDebug(AliLog::kWarning, "Warning: Not implemented");
 208 }
 209
 210 //_____________________________________________________________________________
 211 void AliPerformanceTPC::Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, AliESDfriend *const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend)
 212 {
 213   // Process comparison information
 214   //
 215   if(!esdEvent)
 216   {
 217       AliDebug(AliLog::kError, "esdEvent not available");
 218       return;
 219   }
 220   AliHeader* header = 0;
 221   AliGenEventHeader* genHeader = 0;
 222   AliStack* stack = 0;
 223   TArrayF vtxMC(3);
 224
 225   if(bUseMC)
 226   {
 227     if(!mcEvent) {
 228       AliDebug(AliLog::kError, "mcEvent not available");
 229       return;
 230     }
 231
 232     // get MC event header
 233     header = mcEvent->Header();
 234     if (!header) {
 235       AliDebug(AliLog::kError, "Header not available");
 236       return;
 237     }
 238     // MC particle stack
 239     stack = mcEvent->Stack();
 240     if (!stack) {
 241       AliDebug(AliLog::kError, "Stack not available");
 242       return;
 243     }
 244
 245     // get MC vertex
 246     genHeader = header->GenEventHeader();
 247     if (!genHeader) {
 248       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
 249       return;
 250     }
 251     genHeader->PrimaryVertex(vtxMC);
 252
 253   } // end bUseMC
 254
 255   // use ESD friends
 256   if(bUseESDfriend) {
 257     if(!esdFriend) {
 258       AliDebug(AliLog::kError, "esdFriend not available");
 259       return;
 260     }
 261   }
 262
 263
 264
 265
 266   //  Process events
 267   for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
 268   {
 269     AliESDtrack*track = esdEvent->GetTrack(iTrack);
 270     if(!track) continue;
 271
 272     if(GetAnalysisMode() == 0) ProcessTPC(stack,track);
 273     else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track);
 274     else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track);
 275     else {
 276       printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
 277       return;
 278     }
 279   }
 280 }
 281
 282 //_____________________________________________________________________________
 283 void AliPerformanceTPC::Analyse() {
 284   //
 285   // Analyse comparison information and store output histograms
 286   // in the folder "folderTPC"
 287   //
 288   TH1::AddDirectory(kFALSE);
 289   //TH1F *h=0;
 290   TH2F *h2D=0;
 291   TObjArray *aFolderObj = new TObjArray;
 292
 293   char name[256];
 294   char title[256];
 295   for(Int_t i=0; i<5; i++)
 296   {
 297     for(Int_t j=i+1; j<6; j++)
 298     {
 299       if(j==5) fTPCHisto->GetAxis(5)->SetRangeUser(0.1,10.);
 300       h2D = (TH2F*)fTPCHisto->Projection(i,j);
 301       sprintf(name,"h_tpc_%d_vs_%d",i,j);
 302       h2D->SetName(name);
 303       h2D->GetXaxis()->SetTitle(fTPCHisto->GetAxis(j)->GetTitle());
 304       h2D->GetYaxis()->SetTitle(fTPCHisto->GetAxis(i)->GetTitle());
 305       sprintf(title,"%s vs %s",fTPCHisto->GetAxis(j)->GetTitle(),fTPCHisto->GetAxis(i)->GetTitle());
 306       h2D->SetTitle(title);
 307
 308       if(j==5) h2D->SetBit(TH1::kLogX);
 309       aFolderObj->Add(h2D);
 310     }
 311   }
 312
 313   // export objects to analysis folder
 314   fAnalysisFolder = ExportToFolder(aFolderObj);
 315
 316   // delete only TObjArray
 317   if(aFolderObj) delete aFolderObj;
 318 }
 319
 320 //_____________________________________________________________________________
 321 TFolder* AliPerformanceTPC::ExportToFolder(TObjArray * array)
 322 {
 323   // recreate folder avery time and export objects to new one
 324   //
 325   AliPerformanceTPC * comp=this;
 326   TFolder *folder = comp->GetAnalysisFolder();
 327
 328   TString name, title;
 329   TFolder *newFolder = 0;
 330   Int_t i = 0;
 331   Int_t size = array->GetSize();
 332
 333   if(folder) {
 334      // get name and title from old folder
 335      name = folder->GetName();
 336      title = folder->GetTitle();
 337
 338          // delete old one
 339      delete folder;
 340
 341          // create new one
 342      newFolder = CreateFolder(name.Data(),title.Data());
 343      newFolder->SetOwner();
 344
 345          // add objects to folder
 346      while(i < size) {
 347            newFolder->Add(array->At(i));
 348            i++;
 349          }
 350   }
 351
 352 return newFolder;
 353 }
 354
 355 //_____________________________________________________________________________
 356 Long64_t AliPerformanceTPC::Merge(TCollection* const list)
 357 {
 358   // Merge list of objects (needed by PROOF)
 359
 360   if (!list)
 361   return 0;
 362
 363   if (list->IsEmpty())
 364   return 1;
 365
 366   TIterator* iter = list->MakeIterator();
 367   TObject* obj = 0;
 368
 369   // collection of generated histograms
 370   Int_t count=0;
 371   while((obj = iter->Next()) != 0)
 372   {
 373   AliPerformanceTPC* entry = dynamic_cast<AliPerformanceTPC*>(obj);
 374   if (entry == 0) continue;
 375
 376   fTPCHisto->Add(entry->fTPCHisto);
 377
 378   count++;
 379   }
 380
 381 return count;
 382 }
 383
 384 //_____________________________________________________________________________
 385 TFolder* AliPerformanceTPC::CreateFolder(TString name,TString title) {
 386 // create folder for analysed histograms
 387 //
 388 TFolder *folder = 0;
 389   folder = new TFolder(name.Data(),title.Data());
 390
 391   return folder;
 392 }