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