PWG1/TPC/AliPerformanceDCA.cxx

   1 //------------------------------------------------------------------------------
   2 // Implementation of AliPerformanceDCA 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
   8 // which is a data member of AliPerformanceDCA.
   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   AliPerformanceDCA * compObj = (AliPerformanceDCA*)coutput->FindObject("AliPerformanceDCA");
  20
  21   // Analyse comparison data
  22   compObj->Analyse();
  23
  24   // the output histograms/graphs will be stored in the folder "folderDCA"
  25   compObj->GetAnalysisFolder()->ls("*");
  26
  27   // user can save whole comparison object (or only folder with anlysed histograms)
  28   // in the seperate output file (e.g.)
  29   TFile fout("Analysed_DCA.root","recreate");
  30   compObj->Write(); // compObj->GetAnalysisFolder()->Write();
  31   fout.Close();
  32
  33 */
  34
  35 #include <TAxis.h>
  36 #include <TCanvas.h>
  37 #include <TGraph.h>
  38 #include <TGraph2D.h>
  39 #include <TH1.h>
  40 #include <TH2.h>
  41 #include <TH3.h>
  42 #include <TF1.h>
  43
  44 #include "AliPerformanceDCA.h"
  45 #include "AliESDEvent.h"
  46 #include "AliESDVertex.h"
  47 #include "AliLog.h"
  48 #include "AliMathBase.h"
  49 #include "AliRecInfoCuts.h"
  50 #include "AliMCInfoCuts.h"
  51 #include "AliStack.h"
  52 #include "AliMCEvent.h"
  53 #include "AliTracker.h"
  54 #include "AliHeader.h"
  55 #include "AliGenEventHeader.h"
  56
  57 using namespace std;
  58
  59 ClassImp(AliPerformanceDCA)
  60
  61 //_____________________________________________________________________________
  62 AliPerformanceDCA::AliPerformanceDCA():
  63   AliPerformanceObject("AliPerformanceDCA"),
  64
  65   // DCA histograms
  66   fDCAHisto(0),
  67
  68   // Cuts
  69   fCutsRC(0),
  70   fCutsMC(0),
  71
  72   // histogram folder
  73   fAnalysisFolder(0)
  74 {
  75   // default constructor
  76   Init();
  77 }
  78
  79 //_____________________________________________________________________________
  80 AliPerformanceDCA::AliPerformanceDCA(Char_t* name="AliPerformanceDCA", Char_t* title="AliPerformanceDCA",Int_t analysisMode=0, Bool_t hptGenerator=kFALSE):
  81   AliPerformanceObject(name,title),
  82
  83   // DCA histograms
  84   fDCAHisto(0),
  85
  86   // Cuts
  87   fCutsRC(0),
  88   fCutsMC(0),
  89
  90   // histogram folder
  91   fAnalysisFolder(0)
  92 {
  93   // named constructor
  94
  95   SetAnalysisMode(analysisMode);
  96   SetHptGenerator(hptGenerator);
  97   Init();
  98 }
  99
 100 //_____________________________________________________________________________
 101 AliPerformanceDCA::~AliPerformanceDCA()
 102 {
 103   // destructor
 104   if(fDCAHisto)  delete fDCAHisto; fDCAHisto=0;
 105   if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
 106 }
 107
 108 //_____________________________________________________________________________
 109 void AliPerformanceDCA::Init()
 110 {
 111   // DCA histograms
 112   Int_t nPtBins = 50;
 113   Double_t ptMin = 1.e-2, ptMax = 10.;
 114
 115   Double_t *binsPt = 0;
 116   if (IsHptGenerator())  {
 117     nPtBins = 100; ptMax = 100.;
 118     binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
 119   } else {
 120     binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
 121   }
 122
 123   /*
 124   Int_t nPtBins = 31;
 125    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.};
 126    Double_t ptMin = 0., ptMax = 10.;
 127
 128    if(IsHptGenerator() == kTRUE) {
 129      nPtBins = 100;
 130      ptMin = 0.; ptMax = 100.;
 131    }
 132    */
 133
 134    //dca_r, dca_z, eta, pt
 135    Int_t binsQA[5]    = {100,100,30,nPtBins,144};
 136    Double_t xminQA[5] = {-10.,-10.,-1.5,ptMin,0.};
 137    Double_t xmaxQA[5] = {10.,10.,1.5,ptMax,2*TMath::Pi()};
 138
 139    fDCAHisto = new THnSparseF("fDCAHisto","dca_r:dca_z:eta:pt:phi",5,binsQA,xminQA,xmaxQA);
 140    fDCAHisto->SetBinEdges(3,binsPt);
 141
 142    fDCAHisto->GetAxis(0)->SetTitle("dca_r (cm)");
 143    fDCAHisto->GetAxis(1)->SetTitle("dca_z (cm)");
 144    fDCAHisto->GetAxis(2)->SetTitle("#eta");
 145    fDCAHisto->GetAxis(3)->SetTitle("p_{T} (GeV/c)");
 146    fDCAHisto->GetAxis(4)->SetTitle("phi (rad)");
 147    fDCAHisto->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("folderDCA","Analysis DCA Folder");
 157 }
 158
 159 //_____________________________________________________________________________
 160 void AliPerformanceDCA::ProcessTPC(AliStack* const stack, AliESDtrack *const esdTrack, AliESDEvent* const esdEvent)
 161 {
 162   // Fill DCA comparison information
 163   if(!esdEvent) return;
 164   if(!esdTrack) return;
 165
 166   if( IsUseTrackVertex() )
 167   {
 168     // Relate TPC inner params to prim. vertex
 169     const AliESDVertex *vtxESD = esdEvent->GetPrimaryVertexTracks();
 170     Double_t x[3]; esdTrack->GetXYZ(x);
 171     Double_t b[3]; AliTracker::GetBxByBz(x,b);
 172     Bool_t isOK = esdTrack->RelateToVertexTPCBxByBz(vtxESD, b, kVeryBig);
 173     if(!isOK) return;
 174
 175     /*
 176       // JMT -- recaluclate DCA for HLT if not present
 177       if ( dca[0] == 0. && dca[1] == 0. ) {
 178         track->GetDZ( vtxESD->GetX(), vtxESD->GetY(), vtxESD->GetZ(), esdEvent->GetMagneticField(), dca );
 179       }
 180     */
 181   }
 182
 183   // get TPC inner params at DCA to prim. vertex
 184   const AliExternalTrackParam *track = esdTrack->GetTPCInnerParam();
 185   if(!track) return;
 186
 187   // read from ESD track
 188   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 189   esdTrack->GetImpactParametersTPC(dca,cov);
 190
 191   if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
 192
 193   Double_t vDCAHisto[5]={dca[0],dca[1],track->Eta(),track->Pt(),track->Phi()};
 194   fDCAHisto->Fill(vDCAHisto);
 195
 196   //
 197   // Fill rec vs MC information
 198   //
 199   if(!stack) return;
 200
 201 }
 202
 203 //_____________________________________________________________________________
 204 void AliPerformanceDCA::ProcessTPCITS(AliStack* const stack, AliESDtrack *const esdTrack, AliESDEvent* const esdEvent)
 205 {
 206   // Fill DCA comparison information
 207   if(!esdTrack) return;
 208   if(!esdEvent) return;
 209
 210   if( IsUseTrackVertex() )
 211   {
 212     // Relate TPC inner params to prim. vertex
 213     const AliESDVertex *vtxESD = esdEvent->GetPrimaryVertexTracks();
 214     Double_t x[3]; esdTrack->GetXYZ(x);
 215     Double_t b[3]; AliTracker::GetBxByBz(x,b);
 216     Bool_t isOK = esdTrack->RelateToVertexBxByBz(vtxESD, b, kVeryBig);
 217     if(!isOK) return;
 218
 219     /*
 220       // JMT -- recaluclate DCA for HLT if not present
 221       if ( dca[0] == 0. && dca[1] == 0. ) {
 222         track->GetDZ( vtxESD->GetX(), vtxESD->GetY(), vtxESD->GetZ(), esdEvent->GetMagneticField(), dca );
 223       }
 224     */
 225   }
 226
 227   Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
 228   esdTrack->GetImpactParameters(dca,cov);
 229
 230   if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
 231   if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
 232   if(esdTrack->GetITSclusters(0)<fCutsRC->GetMinNClustersITS()) return;  // min. nb. ITS clusters
 233
 234   Double_t vDCAHisto[5]={dca[0],dca[1],esdTrack->Eta(),esdTrack->Pt(), esdTrack->Phi()};
 235   fDCAHisto->Fill(vDCAHisto);
 236
 237   //
 238   // Fill rec vs MC information
 239   //
 240   if(!stack) return;
 241
 242 }
 243
 244 void AliPerformanceDCA::ProcessConstrained(AliStack* const /*stack*/, AliESDtrack *const /*esdTrack*/)
 245 {
 246   // Fill DCA comparison information
 247
 248   AliDebug(AliLog::kWarning, "Warning: Not implemented");
 249 }
 250
 251 //_____________________________________________________________________________
 252 Long64_t AliPerformanceDCA::Merge(TCollection* const list)
 253 {
 254   // Merge list of objects (needed by PROOF)
 255
 256   if (!list)
 257   return 0;
 258
 259   if (list->IsEmpty())
 260   return 1;
 261
 262   TIterator* iter = list->MakeIterator();
 263   TObject* obj = 0;
 264
 265   // collection of generated histograms
 266   Int_t count=0;
 267   while((obj = iter->Next()) != 0)
 268   {
 269     AliPerformanceDCA* entry = dynamic_cast<AliPerformanceDCA*>(obj);
 270     if (entry == 0) continue;
 271
 272     fDCAHisto->Add(entry->fDCAHisto);
 273     count++;
 274   }
 275
 276 return count;
 277 }
 278
 279 //_____________________________________________________________________________
 280 void AliPerformanceDCA::Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, AliESDfriend *const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend)
 281 {
 282   // Process comparison information
 283   //
 284   if(!esdEvent)
 285   {
 286     Error("Exec","esdEvent not available");
 287     return;
 288   }
 289   AliHeader* header = 0;
 290   AliGenEventHeader* genHeader = 0;
 291   AliStack* stack = 0;
 292   TArrayF vtxMC(3);
 293
 294   if(bUseMC)
 295   {
 296     if(!mcEvent) {
 297       Error("Exec","mcEvent not available");
 298       return;
 299     }
 300     // get MC event header
 301     header = mcEvent->Header();
 302     if (!header) {
 303       Error("Exec","Header not available");
 304       return;
 305     }
 306     // MC particle stack
 307     stack = mcEvent->Stack();
 308     if (!stack) {
 309       Error("Exec","Stack not available");
 310       return;
 311     }
 312     // get MC vertex
 313     genHeader = header->GenEventHeader();
 314     if (!genHeader) {
 315       Error("Exec","Could not retrieve genHeader from Header");
 316       return;
 317     }
 318     genHeader->PrimaryVertex(vtxMC);
 319   }
 320
 321   // use ESD friends
 322   if(bUseESDfriend) {
 323     if(!esdFriend) {
 324       Error("Exec","esdFriend not available");
 325       return;
 326     }
 327   }
 328
 329   // trigger
 330   if(!bUseMC &&GetTriggerClass()) {
 331     Bool_t isEventTriggered = esdEvent->IsTriggerClassFired(GetTriggerClass());
 332     if(!isEventTriggered) return;
 333   }
 334
 335   // get event vertex
 336   const AliESDVertex *vtxESD = NULL;
 337   if( IsUseTrackVertex() )
 338   {
 339     // track vertex
 340     vtxESD = esdEvent->GetPrimaryVertexTracks();
 341   }
 342   else {
 343     // TPC track vertex
 344     vtxESD = esdEvent->GetPrimaryVertexTPC();
 345   }
 346   if(vtxESD && (vtxESD->GetStatus()<=0)) return;
 347
 348   //  Process events
 349   for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
 350   {
 351     AliESDtrack *track = esdEvent->GetTrack(iTrack);
 352     if(!track) continue;
 353
 354     if(GetAnalysisMode() == 0) ProcessTPC(stack,track,esdEvent);
 355     else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track,esdEvent);
 356     else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track);
 357     else {
 358       printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
 359       return;
 360     }
 361   }
 362 }
 363
 364 //_____________________________________________________________________________
 365 void AliPerformanceDCA::Analyse()
 366 {
 367   //
 368   // Analyse comparison information and store output histograms
 369   // in the analysis folder "folderDCA"
 370   //
 371
 372   TH1::AddDirectory(kFALSE);
 373   TH1F *h1D=0;
 374   TH2F *h2D=0;
 375   TObjArray *aFolderObj = new TObjArray;
 376   char title[256];
 377   TObjArray *arr[6] = {0};
 378   TF1 *f1[6] = {0};
 379
 380   // set pt measurable range
 381   //fDCAHisto->GetAxis(3)->SetRangeUser(0.10,10.);
 382
 383   //
 384   h2D = (TH2F*)fDCAHisto->Projection(0,1); // inverse projection convention
 385   h2D->SetName("dca_r_vs_dca_z");
 386   h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(1)->GetTitle());
 387   h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(0)->GetTitle());
 388   sprintf(title,"%s vs %s",fDCAHisto->GetAxis(0)->GetTitle(),fDCAHisto->GetAxis(1)->GetTitle());
 389   h2D->SetTitle(title);
 390   aFolderObj->Add(h2D);
 391
 392   //
 393   h2D = (TH2F*)fDCAHisto->Projection(0,2);
 394   h2D->SetName("dca_r_vs_eta");
 395   h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
 396   h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(0)->GetTitle());
 397   sprintf(title,"%s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(0)->GetTitle());
 398   h2D->SetTitle(title);
 399   aFolderObj->Add(h2D);
 400
 401   //
 402   // mean and rms
 403   //
 404   h1D = MakeStat1D(h2D,0,0);
 405   h1D->SetName("mean_dca_r_vs_eta");
 406   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
 407   h1D->GetYaxis()->SetTitle("mean_dca_r (cm)");
 408   sprintf(title," mean_dca_r (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
 409   h1D->SetTitle(title);
 410   aFolderObj->Add(h1D);
 411
 412   h1D = MakeStat1D(h2D,0,1);
 413   h1D->SetName("rms_dca_r_vs_eta");
 414   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
 415   h1D->GetYaxis()->SetTitle("rms_dca_r (cm)");
 416   sprintf(title," rms_dca_r (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
 417   h1D->SetTitle(title);
 418   aFolderObj->Add(h1D);
 419
 420   //
 421   // fit mean and sigma
 422   //
 423
 424   arr[0] = new TObjArray();
 425   f1[0] = new TF1("gaus","gaus");
 426   h2D->FitSlicesY(f1[0],0,-1,0,"QNR",arr[0]);
 427
 428   h1D = (TH1F*)arr[0]->At(1);
 429   h1D->SetName("fit_mean_dca_r_vs_eta");
 430   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
 431   h1D->GetYaxis()->SetTitle("fit_mean_dca_r (cm)");
 432   sprintf(title," fit_mean_dca_r (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
 433   h1D->SetTitle(title);
 434   aFolderObj->Add(h1D);
 435
 436   h1D = (TH1F*)arr[0]->At(2);
 437   h1D->SetName("res_dca_r_vs_eta");
 438   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
 439   h1D->GetYaxis()->SetTitle("res_dca_r (cm)");
 440   sprintf(title," res_dca_r (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
 441   h1D->SetTitle(title);
 442   aFolderObj->Add(h1D);
 443
 444   //
 445   //
 446   //
 447   h2D = (TH2F*)fDCAHisto->Projection(0,3);
 448   h2D->SetName("dca_r_vs_pt");
 449   h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
 450   h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(0)->GetTitle());
 451   sprintf(title,"%s vs %s",fDCAHisto->GetAxis(0)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
 452   h2D->SetTitle(title);
 453   h2D->SetBit(TH1::kLogX);
 454   aFolderObj->Add(h2D);
 455
 456   h1D = MakeStat1D(h2D,0,0);
 457   h1D->SetName("mean_dca_r_vs_pt");
 458   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
 459   h1D->GetYaxis()->SetTitle("mean_dca_r (cm)");
 460   sprintf(title,"mean_dca_r (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
 461   h1D->SetTitle(title);
 462   h1D->SetBit(TH1::kLogX);
 463   aFolderObj->Add(h1D);
 464
 465   h1D = MakeStat1D(h2D,0,1);
 466   h1D->SetName("rms_dca_r_vs_pt");
 467   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
 468   h1D->GetYaxis()->SetTitle("rms_dca_r (cm)");
 469   sprintf(title,"rms_dca_r (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
 470   h1D->SetTitle(title);
 471   h1D->SetBit(TH1::kLogX);
 472   aFolderObj->Add(h1D);
 473
 474   //
 475   // fit mean and sigma
 476   //
 477
 478   arr[1] = new TObjArray();
 479   f1[1] = new TF1("gaus","gaus");
 480   h2D->FitSlicesY(f1[1],0,-1,0,"QNR",arr[1]);
 481
 482   h1D = (TH1F*)arr[1]->At(1);
 483   h1D->SetName("fit_mean_dca_r_vs_pt");
 484   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
 485   h1D->GetYaxis()->SetTitle("fit_mean_dca_r (cm)");
 486   sprintf(title,"fit_mean_dca_r (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
 487   h1D->SetTitle(title);
 488   h1D->SetBit(TH1::kLogX);
 489   aFolderObj->Add(h1D);
 490
 491   h1D = (TH1F*)arr[1]->At(2);
 492   h1D->SetName("res_dca_r_vs_pt");
 493   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
 494   h1D->GetYaxis()->SetTitle("res_dca_r (cm)");
 495   sprintf(title,"res_dca_r (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
 496   h1D->SetTitle(title);
 497   h1D->SetBit(TH1::kLogX);
 498   aFolderObj->Add(h1D);
 499
 500   //
 501   h2D = (TH2F*)fDCAHisto->Projection(1,2);
 502   h2D->SetName("dca_z_vs_eta");
 503   h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
 504   h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(1)->GetTitle());
 505   sprintf(title,"%s vs %s",fDCAHisto->GetAxis(1)->GetTitle(),fDCAHisto->GetAxis(2)->GetTitle());
 506   h2D->SetTitle(title);
 507   aFolderObj->Add(h2D);
 508
 509   h1D = MakeStat1D(h2D,0,0);
 510   h1D->SetName("mean_dca_z_vs_eta");
 511   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
 512   h1D->GetYaxis()->SetTitle("mean_dca_z (cm)");
 513   sprintf(title,"mean_dca_z (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
 514   h1D->SetTitle(title);
 515   aFolderObj->Add(h1D);
 516
 517   h1D = MakeStat1D(h2D,0,1);
 518   h1D->SetName("rms_dca_z_vs_eta");
 519   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
 520   h1D->GetYaxis()->SetTitle("rms_dca_z (cm)");
 521   sprintf(title,"rms_dca_z (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
 522   h1D->SetTitle(title);
 523   aFolderObj->Add(h1D);
 524
 525   //
 526   // fit mean and sigma
 527   //
 528   arr[2] = new TObjArray();
 529   f1[2] = new TF1("gaus","gaus");
 530   h2D->FitSlicesY(f1[2],0,-1,0,"QNR",arr[2]);
 531
 532   h1D = (TH1F*)arr[2]->At(1);
 533   h1D->SetName("fit_mean_dca_z_vs_eta");
 534   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
 535   h1D->GetYaxis()->SetTitle("fit_mean_dca_z (cm)");
 536   sprintf(title,"fit_mean_dca_z (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
 537   h1D->SetTitle(title);
 538   aFolderObj->Add(h1D);
 539
 540   h1D = (TH1F*)arr[2]->At(2);
 541   h1D->SetName("res_dca_z_vs_eta");
 542   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
 543   h1D->GetYaxis()->SetTitle("res_dca_z (cm)");
 544   sprintf(title,"res_dca_z (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
 545   h1D->SetTitle(title);
 546   aFolderObj->Add(h1D);
 547
 548   //
 549   h2D = (TH2F*)fDCAHisto->Projection(1,3);
 550   h2D->SetName("dca_z_vs_pt");
 551   h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
 552   h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(1)->GetTitle());
 553   sprintf(title,"%s vs %s",fDCAHisto->GetAxis(1)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
 554   h2D->SetTitle(title);
 555   h2D->SetBit(TH1::kLogX);
 556   aFolderObj->Add(h2D);
 557
 558   h1D = MakeStat1D(h2D,0,0);
 559   h1D->SetName("mean_dca_z_vs_pt");
 560   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
 561   h1D->GetYaxis()->SetTitle("mean_dca_z (cm)");
 562   sprintf(title,"mean_dca_z (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
 563   h1D->SetTitle(title);
 564   h1D->SetBit(TH1::kLogX);
 565   aFolderObj->Add(h1D);
 566
 567   h1D = MakeStat1D(h2D,0,1);
 568   h1D->SetName("rms_dca_z_vs_pt");
 569   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
 570   h1D->GetYaxis()->SetTitle("rms_dca_z (cm)");
 571   sprintf(title,"rms_dca_z (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
 572   h1D->SetTitle(title);
 573   h1D->SetBit(TH1::kLogX);
 574   aFolderObj->Add(h1D);
 575
 576   //
 577   // fit mean and sigma
 578   //
 579
 580   arr[3] = new TObjArray();
 581   f1[3] = new TF1("gaus","gaus");
 582   h2D->FitSlicesY(f1[3],0,-1,0,"QNR",arr[3]);
 583
 584   h1D = (TH1F*)arr[3]->At(1);
 585   h1D->SetName("fit_mean_dca_z_vs_pt");
 586   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
 587   h1D->GetYaxis()->SetTitle("fit_mean_dca_z (cm)");
 588   sprintf(title,"fit_mean_dca_z (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
 589   h1D->SetTitle(title);
 590   h1D->SetBit(TH1::kLogX);
 591   aFolderObj->Add(h1D);
 592
 593   h1D = (TH1F*)arr[3]->At(2);
 594   h1D->SetName("res_dca_z_vs_pt");
 595   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
 596   h1D->GetYaxis()->SetTitle("res_dca_z (cm)");
 597   sprintf(title,"res_dca_z (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
 598   h1D->SetTitle(title);
 599   h1D->SetBit(TH1::kLogX);
 600   aFolderObj->Add(h1D);
 601
 602   // A - side
 603   fDCAHisto->GetAxis(2)->SetRangeUser(-1.5,0.0);
 604
 605   h2D = (TH2F*)fDCAHisto->Projection(1,4);
 606   h2D->SetName("dca_z_vs_phi_Aside");
 607   h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(4)->GetTitle());
 608   h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(1)->GetTitle());
 609   sprintf(title,"%s vs %s (A-side)",fDCAHisto->GetAxis(1)->GetTitle(),fDCAHisto->GetAxis(4)->GetTitle());
 610   h2D->SetTitle(title);
 611   aFolderObj->Add(h2D);
 612
 613   h1D = MakeStat1D(h2D,0,0);
 614   h1D->SetName("mean_dca_z_vs_phi_Aside");
 615   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(4)->GetTitle());
 616   h1D->GetYaxis()->SetTitle("mean_dca_z (cm)");
 617   sprintf(title,"mean_dca_z (cm) vs %s (A-side)",fDCAHisto->GetAxis(4)->GetTitle());
 618   h1D->SetTitle(title);
 619   aFolderObj->Add(h1D);
 620
 621   h1D = MakeStat1D(h2D,0,1);
 622   h1D->SetName("rms_dca_z_vs_phi_Aside");
 623   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(4)->GetTitle());
 624   h1D->GetYaxis()->SetTitle("rms_dca_z (cm)");
 625   sprintf(title,"rms_dca_z (cm) vs %s (A-side)",fDCAHisto->GetAxis(4)->GetTitle());
 626   h1D->SetTitle(title);
 627   aFolderObj->Add(h1D);
 628
 629   //
 630   // fit mean and sigma
 631   //
 632   arr[4] = new TObjArray();
 633   f1[4] = new TF1("gaus","gaus");
 634   h2D->FitSlicesY(f1[4],0,-1,0,"QNR",arr[4]);
 635
 636   h1D = (TH1F*)arr[4]->At(1);
 637   h1D->SetName("fit_mean_dca_z_vs_phi_Aside");
 638   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(4)->GetTitle());
 639   h1D->GetYaxis()->SetTitle("fit_mean_dca_z (cm)");
 640   sprintf(title,"fit_mean_dca_z (cm) vs %s (A-side)",fDCAHisto->GetAxis(4)->GetTitle());
 641   h1D->SetTitle(title);
 642   aFolderObj->Add(h1D);
 643
 644   h1D = (TH1F*)arr[4]->At(2);
 645   h1D->SetName("res_dca_z_vs_phi_Aside");
 646   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(4)->GetTitle());
 647   h1D->GetYaxis()->SetTitle("res_dca_z (cm)");
 648   sprintf(title,"res_dca_z (cm) vs %s (A-side)",fDCAHisto->GetAxis(4)->GetTitle());
 649   h1D->SetTitle(title);
 650   aFolderObj->Add(h1D);
 651
 652
 653   // C - side
 654   fDCAHisto->GetAxis(2)->SetRangeUser(0.0,1.5);
 655
 656   h2D = (TH2F*)fDCAHisto->Projection(1,4);
 657   h2D->SetName("dca_z_vs_phi_Cside");
 658   h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(4)->GetTitle());
 659   h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(1)->GetTitle());
 660   sprintf(title,"%s vs %s (C-side)",fDCAHisto->GetAxis(1)->GetTitle(),fDCAHisto->GetAxis(4)->GetTitle());
 661   h2D->SetTitle(title);
 662   aFolderObj->Add(h2D);
 663
 664   h1D = MakeStat1D(h2D,0,0);
 665   h1D->SetName("mean_dca_z_vs_phi_Cside");
 666   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(4)->GetTitle());
 667   h1D->GetYaxis()->SetTitle("mean_dca_z (cm)");
 668   sprintf(title,"mean_dca_z (cm) vs %s (C-side)",fDCAHisto->GetAxis(4)->GetTitle());
 669   h1D->SetTitle(title);
 670   aFolderObj->Add(h1D);
 671
 672   h1D = MakeStat1D(h2D,0,1);
 673   h1D->SetName("rms_dca_z_vs_phi_Cside");
 674   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(4)->GetTitle());
 675   h1D->GetYaxis()->SetTitle("rms_dca_z (cm)");
 676   sprintf(title,"rms_dca_z (cm) vs %s (C-side)",fDCAHisto->GetAxis(4)->GetTitle());
 677   h1D->SetTitle(title);
 678   aFolderObj->Add(h1D);
 679
 680   //
 681   // fit mean and sigma
 682   //
 683   arr[5] = new TObjArray();
 684   f1[5] = new TF1("gaus","gaus");
 685   h2D->FitSlicesY(f1[5],0,-1,0,"QNR",arr[5]);
 686
 687   h1D = (TH1F*)arr[5]->At(1);
 688   h1D->SetName("fit_mean_dca_z_vs_phi_Cside");
 689   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(4)->GetTitle());
 690   h1D->GetYaxis()->SetTitle("fit_mean_dca_z (cm)");
 691   sprintf(title,"fit_mean_dca_z (cm) vs %s (C-side)",fDCAHisto->GetAxis(4)->GetTitle());
 692   h1D->SetTitle(title);
 693   aFolderObj->Add(h1D);
 694
 695   h1D = (TH1F*)arr[5]->At(2);
 696   h1D->SetName("res_dca_z_vs_phi_Cside");
 697   h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(4)->GetTitle());
 698   h1D->GetYaxis()->SetTitle("res_dca_z (cm)");
 699   sprintf(title,"res_dca_z (cm) vs %s (C-side)",fDCAHisto->GetAxis(4)->GetTitle());
 700   h1D->SetTitle(title);
 701   aFolderObj->Add(h1D);
 702
 703   // export objects to analysis folder
 704   fAnalysisFolder = ExportToFolder(aFolderObj);
 705
 706   // delete only TObjArray
 707   if(aFolderObj) delete aFolderObj;
 708 }
 709
 710 //_____________________________________________________________________________
 711 TH1F* AliPerformanceDCA::MakeStat1D(TH2 *hist, Int_t delta0, Int_t type)
 712 {
 713   // Return TH1F histogram
 714   // delta - number of bins to integrate
 715   // with mean (type == 0) or RMS (type==1)
 716
 717   char hname[256];
 718   const char* suffix = "_stat1d";
 719   sprintf(hname,"%s%s",hist->GetName(),suffix);
 720   TAxis* xaxis = hist->GetXaxis();
 721   Int_t  nbinx = xaxis->GetNbins();
 722
 723   TH1F *hnew = (TH1F*)hist->ProjectionX()->Clone();
 724   hnew->SetName(hname);
 725
 726   char name[256];
 727   for (Int_t ix=0; ix<=nbinx;ix++) {
 728     sprintf(name,"%s_%d",hist->GetName(),ix);
 729     TH1 *projection = hist->ProjectionY(name,ix-delta0,ix+delta0);
 730
 731     Float_t stat= 0., stat_err =0.;
 732     if (type==0) { stat = projection->GetMean(); stat_err = projection->GetMeanError(); }
 733     if (type==1) { stat = projection->GetRMS(); stat_err = projection->GetRMSError(); }
 734
 735     hnew->SetBinContent(ix, stat);
 736     hnew->SetBinError(ix, stat_err);
 737   }
 738
 739 return hnew;
 740 }
 741
 742 //_____________________________________________________________________________
 743 TH2F* AliPerformanceDCA::MakeStat2D(TH3 *hist, Int_t delta0, Int_t delta1, Int_t type)
 744 {
 745   // Return TH1F histogram
 746   // delta0 - number of bins to integrate in x
 747   // delta1 - number of bins to integrate in y
 748   // with mean (type==0) or RMS (type==1)
 749
 750   char hname[256];
 751   const char* suffix = "_stat2d";
 752   sprintf(hname,"%s%s",hist->GetName(),suffix);
 753
 754   TAxis* xaxis = hist->GetXaxis();
 755   Int_t  nbinx = xaxis->GetNbins();
 756
 757   TH2F *hnew = (TH2F*)hist->Project3D("yx")->Clone();
 758   hnew->SetName(hname);
 759
 760   TAxis* yaxis = hist->GetYaxis();
 761   Int_t  nbiny = yaxis->GetNbins();
 762
 763   char name[256];
 764   for (Int_t ix=0; ix<=nbinx;ix++) {
 765     for (Int_t iy=0; iy<=nbiny;iy++) {
 766       sprintf(name,"%s_%d_%d",hist->GetName(),ix,iy);
 767       TH1 *projection = hist->ProjectionZ(name,ix-delta0,ix+delta0,iy-delta1,iy+delta1);
 768
 769       Float_t stat= 0., stat_err =0.;
 770       if (type==0) { stat = projection->GetMean(); stat_err = projection->GetMeanError(); }
 771       if (type==1) { stat = projection->GetRMS(); stat_err = projection->GetRMSError(); }
 772
 773       hnew->SetBinContent(ix,iy,stat);
 774       hnew->SetBinError(ix,iy,stat_err);
 775     }
 776   }
 777
 778 return hnew;
 779 }
 780
 781 //_____________________________________________________________________________
 782 TFolder* AliPerformanceDCA::ExportToFolder(TObjArray * array)
 783 {
 784   // recreate folder avery time and export objects to new one
 785   //
 786   AliPerformanceDCA * comp=this;
 787   TFolder *folder = comp->GetAnalysisFolder();
 788
 789   TString name, title;
 790   TFolder *newFolder = 0;
 791   Int_t i = 0;
 792   Int_t size = array->GetSize();
 793
 794   if(folder) {
 795      // get name and title from old folder
 796      name = folder->GetName();
 797      title = folder->GetTitle();
 798
 799          // delete old one
 800      delete folder;
 801
 802          // create new one
 803      newFolder = CreateFolder(name.Data(),title.Data());
 804      newFolder->SetOwner();
 805
 806          // add objects to folder
 807      while(i < size) {
 808            newFolder->Add(array->At(i));
 809            i++;
 810          }
 811   }
 812
 813 return newFolder;
 814 }
 815
 816 //_____________________________________________________________________________
 817 TFolder* AliPerformanceDCA::CreateFolder(TString name,TString title) {
 818 // create folder for analysed histograms
 819 TFolder *folder = 0;
 820   folder = new TFolder(name.Data(),title.Data());
 821
 822   return folder;
 823 }