PWG1/TPC/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,90,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     Error("Exec","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       Error("Exec","mcEvent not available");
 229       return;
 230     }
 231     // get MC event header
 232     header = mcEvent->Header();
 233     if (!header) {
 234       Error("Exec","Header not available");
 235       return;
 236     }
 237     // MC particle stack
 238     stack = mcEvent->Stack();
 239     if (!stack) {
 240       Error("Exec","Stack not available");
 241       return;
 242     }
 243     // get MC vertex
 244     genHeader = header->GenEventHeader();
 245     if (!genHeader) {
 246       Error("Exec","Could not retrieve genHeader from Header");
 247       return;
 248     }
 249     genHeader->PrimaryVertex(vtxMC);
 250   }
 251
 252   // use ESD friends
 253   if(bUseESDfriend) {
 254     if(!esdFriend) {
 255       Error("Exec","esdFriend not available");
 256       return;
 257     }
 258   }
 259
 260   //  Process events
 261   for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
 262   {
 263     AliESDtrack*track = esdEvent->GetTrack(iTrack);
 264     if(!track) continue;
 265
 266     if(GetAnalysisMode() == 0) ProcessTPC(stack,track);
 267     else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track);
 268     else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track);
 269     else {
 270       printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
 271       return;
 272     }
 273   }
 274 }
 275
 276 //_____________________________________________________________________________
 277 void AliPerformanceTPC::Analyse() {
 278   //
 279   // Analyse comparison information and store output histograms
 280   // in the folder "folderTPC"
 281   //
 282   TH1::AddDirectory(kFALSE);
 283   //TH1F *h=0;
 284   TH2F *h2D=0;
 285   TObjArray *aFolderObj = new TObjArray;
 286
 287   char name[256];
 288   char title[256];
 289   for(Int_t i=0; i<5; i++)
 290   {
 291     for(Int_t j=i+1; j<6; j++)
 292     {
 293       if(j==5) fTPCHisto->GetAxis(5)->SetRangeUser(0.1,10.);
 294       h2D = (TH2F*)fTPCHisto->Projection(i,j);
 295       sprintf(name,"h_tpc_%d_vs_%d",i,j);
 296       h2D->SetName(name);
 297       h2D->GetXaxis()->SetTitle(fTPCHisto->GetAxis(j)->GetTitle());
 298       h2D->GetYaxis()->SetTitle(fTPCHisto->GetAxis(i)->GetTitle());
 299       sprintf(title,"%s vs %s",fTPCHisto->GetAxis(j)->GetTitle(),fTPCHisto->GetAxis(i)->GetTitle());
 300       h2D->SetTitle(title);
 301
 302       if(j==5) h2D->SetBit(TH1::kLogX);
 303       aFolderObj->Add(h2D);
 304     }
 305   }
 306
 307   // export objects to analysis folder
 308   fAnalysisFolder = ExportToFolder(aFolderObj);
 309
 310   // delete only TObjArray
 311   if(aFolderObj) delete aFolderObj;
 312 }
 313
 314 //_____________________________________________________________________________
 315 TFolder* AliPerformanceTPC::ExportToFolder(TObjArray * array)
 316 {
 317   // recreate folder avery time and export objects to new one
 318   //
 319   AliPerformanceTPC * comp=this;
 320   TFolder *folder = comp->GetAnalysisFolder();
 321
 322   TString name, title;
 323   TFolder *newFolder = 0;
 324   Int_t i = 0;
 325   Int_t size = array->GetSize();
 326
 327   if(folder) {
 328      // get name and title from old folder
 329      name = folder->GetName();
 330      title = folder->GetTitle();
 331
 332          // delete old one
 333      delete folder;
 334
 335          // create new one
 336      newFolder = CreateFolder(name.Data(),title.Data());
 337      newFolder->SetOwner();
 338
 339          // add objects to folder
 340      while(i < size) {
 341            newFolder->Add(array->At(i));
 342            i++;
 343          }
 344   }
 345
 346 return newFolder;
 347 }
 348
 349 //_____________________________________________________________________________
 350 Long64_t AliPerformanceTPC::Merge(TCollection* const list)
 351 {
 352   // Merge list of objects (needed by PROOF)
 353
 354   if (!list)
 355   return 0;
 356
 357   if (list->IsEmpty())
 358   return 1;
 359
 360   TIterator* iter = list->MakeIterator();
 361   TObject* obj = 0;
 362
 363   // collection of generated histograms
 364   Int_t count=0;
 365   while((obj = iter->Next()) != 0)
 366   {
 367   AliPerformanceTPC* entry = dynamic_cast<AliPerformanceTPC*>(obj);
 368   if (entry == 0) continue;
 369
 370   fTPCHisto->Add(entry->fTPCHisto);
 371
 372   count++;
 373   }
 374
 375 return count;
 376 }
 377
 378 //_____________________________________________________________________________
 379 TFolder* AliPerformanceTPC::CreateFolder(TString name,TString title) {
 380 // create folder for analysed histograms
 381 //
 382 TFolder *folder = 0;
 383   folder = new TFolder(name.Data(),title.Data());
 384
 385   return folder;
 386 }