PWG1/TPC/AliPerformanceDEdx.cxx

   1 //------------------------------------------------------------------------------
   2 // Implementation of AliPerformanceDEdx 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 of AliPerformanceDEdx.
   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   TFile f("Output.root");
  19   //AliPerformanceDEdx * compObj = (AliPerformanceDEdx*)f.Get("AliPerformanceDEdx");
  20   AliPerformanceDEdx * compObj = (AliPerformanceDEdx*)coutput->FindObject("AliPerformanceDEdx");
  21
  22   // Analyse comparison data
  23   compObj->Analyse();
  24
  25   // the output histograms/graphs will be stored in the folder "folderDEdx"
  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_DEdx.root"."recreate");
  31   compObj->Write(); // compObj->GetAnalysisFolder()->Write();
  32   fout.Close();
  33
  34 */
  35
  36 #include <TDirectory.h>
  37 #include <TAxis.h>
  38 #include <TCanvas.h>
  39 #include <TH1.h>
  40 #include <TH2.h>
  41 #include <TF1.h>
  42
  43 #include "AliPerformanceDEdx.h"
  44 #include "AliESDEvent.h"
  45 #include "AliMCEvent.h"
  46 #include "AliESDtrack.h"
  47 #include "AliStack.h"
  48 #include "AliLog.h"
  49 #include "AliMCInfoCuts.h"
  50 #include "AliMathBase.h"
  51 #include "AliRecInfoCuts.h"
  52 #include "AliTreeDraw.h"
  53 #include "AliHeader.h"
  54 #include "AliGenEventHeader.h"
  55
  56 using namespace std;
  57
  58 ClassImp(AliPerformanceDEdx)
  59
  60 //_____________________________________________________________________________
  61 AliPerformanceDEdx::AliPerformanceDEdx():
  62   AliPerformanceObject("AliPerformanceDEdx"),
  63
  64   // dEdx
  65   fDeDxHisto(0),
  66
  67   // Cuts
  68   fCutsRC(0),
  69   fCutsMC(0),
  70
  71   // histogram folder
  72   fAnalysisFolder(0)
  73 {
  74   // default constructor
  75   Init();
  76 }
  77
  78 //_____________________________________________________________________________
  79 AliPerformanceDEdx::AliPerformanceDEdx(Char_t* name="AliPerformanceDEdx", Char_t* title="AliPerformanceDEdx",Int_t analysisMode=0, Bool_t hptGenerator=kFALSE):
  80   AliPerformanceObject(name,title),
  81
  82   // dEdx
  83   fDeDxHisto(0),
  84
  85   // Cuts
  86   fCutsRC(0),
  87   fCutsMC(0),
  88
  89   // histogram folder
  90   fAnalysisFolder(0)
  91 {
  92   // named constructor
  93
  94   SetAnalysisMode(analysisMode);
  95   SetHptGenerator(hptGenerator);
  96   Init();
  97 }
  98
  99
 100 //_____________________________________________________________________________
 101 AliPerformanceDEdx::~AliPerformanceDEdx()
 102 {
 103   // destructor
 104   if(fDeDxHisto)  delete fDeDxHisto; fDeDxHisto=0;
 105   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
 106 }
 107
 108 //_____________________________________________________________________________
 109 void AliPerformanceDEdx::Init()
 110 {
 111   // Init histograms
 112
 113   // TPC dEdx
 114   // set pt bins
 115   Int_t nPBins = 50;
 116   Double_t pMin = 1.e-2, pMax = 10.;
 117
 118   Double_t *binsP = 0;
 119   if (IsHptGenerator())  {
 120     nPBins = 100; pMax = 100.;
 121     binsP = CreateLogAxis(nPBins,pMin,pMax);
 122   } else {
 123     binsP = CreateLogAxis(nPBins,pMin,pMax);
 124   }
 125
 126
 127   /*
 128   Int_t nPBins = 31;
 129     Double_t binsP[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.};
 130     Double_t pMin = 0., pMax = 10.;
 131
 132     if(IsHptGenerator() == kTRUE) {
 133       nPBins = 100;
 134       pMin = 0.; pMax = 100.;
 135     }
 136    */
 137
 138    //dedx:alpha:y:z:snp:tgl:ncls:p
 139    //dedx:phi:y:z:snp:tgl:ncls:p
 140    //Int_t binsQA[8]    = {300, 50, 50,  50, 50, 50, 80, nPBins};
 141    //Double_t xminQA[8] = {0, -4,-20,-250, -1, -2, 0, pMin};
 142    //Double_t xmaxQA[8] = {300, 4, 20, 250,  1,  2, 160, pMax};
 143    Int_t binsQA[8]    = {300, 144, 50,  50, 50, 50, 80, nPBins};
 144    Double_t xminQA[8] = {0, -TMath::Pi(),-20,-250, -1, -2, 0, pMin};
 145    Double_t xmaxQA[8] = {300, TMath::Pi(), 20, 250,  1,  2, 160, pMax};
 146
 147    //fDeDxHisto = new THnSparseF("fDeDxHisto","dedx:alpha:y:z:snp:tgl:ncls:momentum",8,binsQA,xminQA,xmaxQA);
 148    fDeDxHisto = new THnSparseF("fDeDxHisto","dedx:phi:y:z:snp:tgl:ncls:momentum",8,binsQA,xminQA,xmaxQA);
 149    fDeDxHisto->SetBinEdges(7,binsP);
 150
 151    fDeDxHisto->GetAxis(0)->SetTitle("dedx (a.u.)");
 152    fDeDxHisto->GetAxis(1)->SetTitle("#phi (rad)");
 153    fDeDxHisto->GetAxis(2)->SetTitle("y (cm)");
 154    fDeDxHisto->GetAxis(3)->SetTitle("z (cm)");
 155    fDeDxHisto->GetAxis(4)->SetTitle("sin#phi");
 156    fDeDxHisto->GetAxis(5)->SetTitle("tan#lambda");
 157    fDeDxHisto->GetAxis(6)->SetTitle("ncls");
 158    fDeDxHisto->GetAxis(7)->SetTitle("p (GeV/c)");
 159    fDeDxHisto->Sumw2();
 160
 161    // Init cuts
 162    if(!fCutsMC)
 163      AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
 164    if(!fCutsRC)
 165      AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
 166
 167    // init folder
 168    fAnalysisFolder = CreateFolder("folderDEdx","Analysis de/dx Folder");
 169 }
 170
 171 //_____________________________________________________________________________
 172 void AliPerformanceDEdx::ProcessTPC(AliStack* const /*stack*/, AliESDtrack *const /*esdTrack*/)
 173 {
 174   // Fill dE/dx  comparison information
 175   AliDebug(AliLog::kWarning, "Warning: Not implemented");
 176 }
 177
 178 //_____________________________________________________________________________
 179 void AliPerformanceDEdx::ProcessInnerTPC(AliStack* const stack, AliESDtrack *const esdTrack)
 180 {
 181   if(!esdTrack) return;
 182
 183   const AliExternalTrackParam *innerParam =  esdTrack->GetInnerParam();
 184   if(!innerParam) return;
 185
 186   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 187   esdTrack->GetImpactParametersTPC(dca,cov);
 188
 189   if((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
 190
 191   //
 192   // select primaries
 193   //
 194   Double_t dcaToVertex = -1;
 195   if( fCutsRC->GetDCAToVertex2D() )
 196   {
 197       dcaToVertex = TMath::Sqrt(dca[0]*dca[0]/fCutsRC->GetMaxDCAToVertexXY()/fCutsRC->GetMaxDCAToVertexXY()                    + dca[1]*dca[1]/fCutsRC->GetMaxDCAToVertexZ()/fCutsRC->GetMaxDCAToVertexZ());
 198   }
 199   if(fCutsRC->GetDCAToVertex2D() && dcaToVertex > 1) return;
 200   if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[0]) > fCutsRC->GetMaxDCAToVertexXY()) return;
 201   if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[1]) > fCutsRC->GetMaxDCAToVertexZ()) return;
 202
 203   Float_t dedx = esdTrack->GetTPCsignal();
 204   Int_t ncls = esdTrack->GetTPCNcls();
 205
 206   Double_t pt = innerParam->Pt();
 207   Double_t lam = TMath::ATan2(innerParam->Pz(),innerParam->Pt());
 208   Double_t p = pt/TMath::Cos(lam);
 209   //Double_t alpha = innerParam->GetAlpha();
 210   Double_t phi = TMath::ATan2(innerParam->Py(),innerParam->Px());
 211   //if(phi<0.) phi += 2.*TMath::Phi();
 212   Double_t y = innerParam->GetY();
 213   Double_t z = innerParam->GetZ();
 214   Double_t snp = innerParam->GetSnp();
 215   Double_t tgl = innerParam->GetTgl();
 216
 217   //Double_t vDeDxHisto[8] = {dedx,alpha,y,z,snp,tgl,ncls,p};
 218   Double_t vDeDxHisto[8] = {dedx,phi,y,z,snp,tgl,ncls,p};
 219   fDeDxHisto->Fill(vDeDxHisto);
 220
 221   if(!stack) return;
 222 }
 223
 224 //_____________________________________________________________________________
 225 void AliPerformanceDEdx::ProcessTPCITS(AliStack* const /*stack*/, AliESDtrack *const /*esdTrack*/)
 226 {
 227   // Fill dE/dx  comparison information
 228
 229    AliDebug(AliLog::kWarning, "Warning: Not implemented");
 230 }
 231
 232 //_____________________________________________________________________________
 233 void AliPerformanceDEdx::ProcessConstrained(AliStack* const /*stack*/, AliESDtrack *const /*esdTrack*/)
 234 {
 235   // Fill dE/dx  comparison information
 236
 237    AliDebug(AliLog::kWarning, "Warning: Not implemented");
 238 }
 239
 240 //_____________________________________________________________________________
 241 Long64_t AliPerformanceDEdx::Merge(TCollection* const list)
 242 {
 243   // Merge list of objects (needed by PROOF)
 244
 245   if (!list)
 246   return 0;
 247
 248   if (list->IsEmpty())
 249   return 1;
 250
 251   TIterator* iter = list->MakeIterator();
 252   TObject* obj = 0;
 253
 254   // collection of generated histograms
 255   Int_t count=0;
 256   while((obj = iter->Next()) != 0)
 257   {
 258     AliPerformanceDEdx* entry = dynamic_cast<AliPerformanceDEdx*>(obj);
 259     if (entry == 0) continue;
 260
 261     fDeDxHisto->Add(entry->fDeDxHisto);
 262     count++;
 263   }
 264
 265 return count;
 266 }
 267
 268 //_____________________________________________________________________________
 269 void AliPerformanceDEdx::Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, AliESDfriend *const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend)
 270 {
 271   // Process comparison information
 272   //
 273   if(!esdEvent)
 274   {
 275       AliDebug(AliLog::kError, "esdEvent not available");
 276       return;
 277   }
 278   AliHeader* header = 0;
 279   AliGenEventHeader* genHeader = 0;
 280   AliStack* stack = 0;
 281   TArrayF vtxMC(3);
 282
 283   if(bUseMC)
 284   {
 285     if(!mcEvent) {
 286       AliDebug(AliLog::kError, "mcEvent not available");
 287       return;
 288     }
 289
 290     // get MC event header
 291     header = mcEvent->Header();
 292     if (!header) {
 293       AliDebug(AliLog::kError, "Header not available");
 294       return;
 295     }
 296     // MC particle stack
 297     stack = mcEvent->Stack();
 298     if (!stack) {
 299       AliDebug(AliLog::kError, "Stack not available");
 300       return;
 301     }
 302
 303     // get MC vertex
 304     genHeader = header->GenEventHeader();
 305     if (!genHeader) {
 306       AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
 307       return;
 308     }
 309     genHeader->PrimaryVertex(vtxMC);
 310
 311   } // end bUseMC
 312
 313   // use ESD friends
 314   if(bUseESDfriend) {
 315     if(!esdFriend) {
 316       AliDebug(AliLog::kError, "esdFriend not available");
 317       return;
 318     }
 319   }
 320
 321   // trigger
 322   if(!bUseMC) {
 323     Bool_t isEventTriggered = esdEvent->IsTriggerClassFired(GetTriggerClass());
 324     if(!isEventTriggered) return;
 325   }
 326
 327   // get TPC event vertex
 328   const AliESDVertex *vtxESD = esdEvent->GetPrimaryVertexTPC();
 329   if(vtxESD && (vtxESD->GetStatus()<=0)) return;
 330
 331   //  Process events
 332   for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
 333   {
 334     AliESDtrack *track = esdEvent->GetTrack(iTrack);
 335     if(!track) continue;
 336
 337     if(GetAnalysisMode() == 0) ProcessTPC(stack,track);
 338     else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track);
 339     else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track);
 340     else if(GetAnalysisMode() == 3) ProcessInnerTPC(stack,track);
 341     else {
 342       printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
 343       return;
 344     }
 345   }
 346 }
 347
 348 //_____________________________________________________________________________
 349 void AliPerformanceDEdx::Analyse()
 350 {
 351   // Analyze comparison information and store output histograms
 352   // in the folder "folderDEdx"
 353   //
 354   TH1::AddDirectory(kFALSE);
 355   TH1F *h1D=0;
 356   TH2F *h2D=0;
 357   TObjArray *aFolderObj = new TObjArray;
 358
 359   char name[256];
 360   char title[256];
 361
 362   for(Int_t i=1; i<8; i++) {
 363     //
 364     h2D = (TH2F*)fDeDxHisto->Projection(0,i);
 365
 366     sprintf(name,"h_dedx_%d_vs_%d",0,i);
 367     h2D->SetName(name);
 368     sprintf(title,"%s vs %s",fDeDxHisto->GetAxis(0)->GetTitle(),fDeDxHisto->GetAxis(i)->GetTitle());
 369     h2D->SetTitle(title);
 370     h2D->GetXaxis()->SetTitle(fDeDxHisto->GetAxis(i)->GetTitle());
 371     h2D->GetYaxis()->SetTitle(fDeDxHisto->GetAxis(0)->GetTitle());
 372
 373     if(i==7) h2D->SetBit(TH1::kLogX);
 374     aFolderObj->Add(h2D);
 375   }
 376
 377   // resolution histograms for mips
 378   //dedx:phi:y:z:snp:tgl:ncls:p
 379   fDeDxHisto->GetAxis(2)->SetRangeUser(-15.,14.999);
 380   fDeDxHisto->GetAxis(3)->SetRangeUser(-120.,119.999);
 381   fDeDxHisto->GetAxis(4)->SetRangeUser(-0.4, 0.399);
 382   fDeDxHisto->GetAxis(5)->SetRangeUser(-0.9,0.89);
 383   fDeDxHisto->GetAxis(6)->SetRangeUser(60.,160.);
 384   fDeDxHisto->GetAxis(7)->SetRangeUser(0.4,0.499);
 385
 386   h1D=(TH1F*)fDeDxHisto->Projection(0);
 387   h1D->SetName("dedx_mips");
 388
 389   h1D->SetTitle("dedx_mips");
 390   h1D->GetXaxis()->SetTitle(fDeDxHisto->GetAxis(0)->GetTitle());
 391   aFolderObj->Add(h1D);
 392
 393   //
 394   TObjArray *arr[7] = {0};
 395   TF1 *f1[7] = {0};
 396
 397   for(Int_t i=1; i<8; i++)
 398   {
 399     arr[i] = new TObjArray;
 400     f1[i] = new TF1("gaus","gaus");
 401     //printf("i %d \n",i);
 402
 403     h2D = (TH2F*)fDeDxHisto->Projection(0,i);
 404
 405     f1[i]->SetRange(40,60); // should be pion peak
 406     h2D->FitSlicesY(f1[i],0,-1,10,"QNR",arr[i]); // gaus fit of pion peak
 407
 408     h1D = (TH1F*)arr[i]->At(1);
 409     sprintf(name,"mean_dedx_mips_vs_%d",i);
 410     h1D->SetName(name);
 411     sprintf(title,"%s vs %s","mean_dedx_mips (a.u.)",fDeDxHisto->GetAxis(i)->GetTitle());
 412     h1D->SetTitle(title);
 413     h1D->GetXaxis()->SetTitle(fDeDxHisto->GetAxis(i)->GetTitle());
 414     h1D->GetYaxis()->SetTitle("mean_dedx_mips (a.u.)");
 415     //h1D->SetMinimum(40);
 416     //h1D->SetMaximum(60);
 417
 418     aFolderObj->Add(h1D);
 419
 420     h1D = (TH1F*)arr[i]->At(2);
 421     sprintf(name,"res_dedx_mips_vs_%d",i);
 422     h1D->SetName(name);
 423     sprintf(title,"%s vs %s","res_dedx_mips (a.u)",fDeDxHisto->GetAxis(i)->GetTitle());
 424     h1D->SetTitle(title);
 425     h1D->GetXaxis()->SetTitle(fDeDxHisto->GetAxis(i)->GetTitle());
 426     h1D->GetYaxis()->SetTitle("res_dedx_mips (a.u.)");
 427     //h1D->SetMinimum(0);
 428     //h1D->SetMaximum(6);
 429
 430     aFolderObj->Add(h1D);
 431   }
 432
 433   // export objects to analysis folder
 434   fAnalysisFolder = ExportToFolder(aFolderObj);
 435
 436   // delete only TObjArray
 437   if(aFolderObj) delete aFolderObj;
 438 }
 439
 440 //_____________________________________________________________________________
 441 TFolder* AliPerformanceDEdx::ExportToFolder(TObjArray * array)
 442 {
 443   // recreate folder avery time and export objects to new one
 444   //
 445   AliPerformanceDEdx * comp=this;
 446   TFolder *folder = comp->GetAnalysisFolder();
 447
 448   TString name, title;
 449   TFolder *newFolder = 0;
 450   Int_t i = 0;
 451   Int_t size = array->GetSize();
 452
 453   if(folder) {
 454      // get name and title from old folder
 455      name = folder->GetName();
 456      title = folder->GetTitle();
 457
 458          // delete old one
 459      delete folder;
 460
 461          // create new one
 462      newFolder = CreateFolder(name.Data(),title.Data());
 463      newFolder->SetOwner();
 464
 465          // add objects to folder
 466      while(i < size) {
 467            newFolder->Add(array->At(i));
 468            i++;
 469          }
 470   }
 471
 472 return newFolder;
 473 }
 474
 475
 476 //_____________________________________________________________________________
 477 TFolder* AliPerformanceDEdx::CreateFolder(TString name,TString title) {
 478 // create folder for analysed histograms
 479 TFolder *folder = 0;
 480   folder = new TFolder(name.Data(),title.Data());
 481
 482   return folder;
 483 }
 484