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