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