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