[u/mrichter/AliRoot.git] / PWG1 / AliPerformanceDCA.cxx

//------------------------------------------------------------------------------
// Implementation of AliPerformanceDCA class. It keeps information from 
// comparison of reconstructed and MC particle tracks. In addtion, 
// it keeps selection cuts used during comparison. The comparison 
// information is stored in the ROOT histograms. Analysis of these 
// histograms can be done by using Analyse() class function. The result of 
// the analysis (histograms/graphs) are stored in the folder
// which is a data member of AliPerformanceDCA.
//  
// Author: J.Otwinowski 04/02/2008 
//------------------------------------------------------------------------------

/*
 
  // after running comparison task, read the file, and get component
  gROOT->LoadMacro("$ALICE_ROOT/PWG1/Macros/LoadMyLibs.C");
  LoadMyLibs();
  TFile f("Output.root");
  AliPerformanceDCA * compObj = (AliPerformanceDCA*)coutput->FindObject("AliPerformanceDCA");

  // Analyse comparison data
  compObj->Analyse();

  // the output histograms/graphs will be stored in the folder "folderDCA" 
  compObj->GetAnalysisFolder()->ls("*");
 
  // user can save whole comparison object (or only folder with anlysed histograms) 
  // in the seperate output file (e.g.)
  TFile fout("Analysed_DCA.root","recreate");
  compObj->Write(); // compObj->GetAnalysisFolder()->Write();
  fout.Close();

*/

#include <TAxis.h>
#include <TCanvas.h>
#include <TGraph.h>
#include <TGraph2D.h>
#include <TH1.h>
#include <TH2.h>
#include <TH3.h>

#include "AliPerformanceDCA.h" 
#include "AliESDEvent.h"   
#include "AliESDVertex.h" 
#include "AliLog.h" 
#include "AliMathBase.h"
#include "AliRecInfoCuts.h" 
#include "AliMCInfoCuts.h" 
#include "AliStack.h" 
#include "AliMCEvent.h" 
#include "AliTracker.h"   
#include "AliHeader.h"   
#include "AliGenEventHeader.h"   

using namespace std;

ClassImp(AliPerformanceDCA)

//_____________________________________________________________________________
AliPerformanceDCA::AliPerformanceDCA():
  AliPerformanceObject("AliPerformanceDCA"),

  // DCA histograms
  fDCAHisto(0),

  // Cuts 
  fCutsRC(0), 
  fCutsMC(0),  

  // histogram folder 
  fAnalysisFolder(0)
{
  // default constructor	
  Init();
}

//_____________________________________________________________________________
AliPerformanceDCA::AliPerformanceDCA(Char_t* name="AliPerformanceDCA", Char_t* title="AliPerformanceDCA",Int_t analysisMode=0, Bool_t hptGenerator=kFALSE):
  AliPerformanceObject(name,title),

  // DCA histograms
  fDCAHisto(0),

  // Cuts 
  fCutsRC(0), 
  fCutsMC(0),  

  // histogram folder 
  fAnalysisFolder(0)
{
  // named constructor	 

  SetAnalysisMode(analysisMode);
  SetHptGenerator(hptGenerator);
  Init();
}

//_____________________________________________________________________________
AliPerformanceDCA::~AliPerformanceDCA()
{
  // destructor
  if(fDCAHisto)  delete fDCAHisto; fDCAHisto=0; 
  if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
}

//_____________________________________________________________________________
void AliPerformanceDCA::Init()
{
  // DCA histograms
  Int_t nPtBins = 50;
  Double_t ptMin = 1.e-2, ptMax = 10.;

  Double_t *binsPt = 0;
  if (IsHptGenerator())  { 
    nPtBins = 100; ptMax = 100.;
    binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
  } else {
    binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
  }

  /*
  Int_t nPtBins = 31;
   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.};
   Double_t ptMin = 0., ptMax = 10.;

   if(IsHptGenerator() == kTRUE) {
     nPtBins = 100;
     ptMin = 0.; ptMax = 100.;
   }
   */

   //dca_r, dca_z, eta, pt
   Int_t binsQA[4]    = {100,100,30,nPtBins};
   Double_t xminQA[4] = {-10.,-10.,-1.5,ptMin};
   Double_t xmaxQA[4] = {10.,10.,1.5,ptMax};

   fDCAHisto = new THnSparseF("fDCAHisto","dca_r:dca_z:eta:pt",4,binsQA,xminQA,xmaxQA);
   fDCAHisto->SetBinEdges(3,binsPt);

   fDCAHisto->GetAxis(0)->SetTitle("dca_r (cm)");
   fDCAHisto->GetAxis(1)->SetTitle("dca_z (cm)");
   fDCAHisto->GetAxis(2)->SetTitle("#eta");
   fDCAHisto->GetAxis(3)->SetTitle("p_{T} (GeV/c)");
   fDCAHisto->Sumw2();

  // init cuts
  if(!fCutsMC) 
    AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
  if(!fCutsRC) 
    AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
 
  // init folder
  fAnalysisFolder = CreateFolder("folderDCA","Analysis DCA Folder");
}

//_____________________________________________________________________________
void AliPerformanceDCA::ProcessTPC(AliStack* const stack, AliESDtrack *const esdTrack)
{
  // Fill DCA comparison information
  if(!esdTrack) return;

  const AliExternalTrackParam *track = esdTrack->GetTPCInnerParam();
  if(!track) return;

  Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
  esdTrack->GetImpactParametersTPC(dca,cov);

  if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters  
 
  Double_t vDCAHisto[4]={dca[0],dca[1],track->Eta(),track->Pt()};
  fDCAHisto->Fill(vDCAHisto);

  //
  // Fill rec vs MC information
  //
  if(!stack) return;

 }

//_____________________________________________________________________________
void AliPerformanceDCA::ProcessTPCITS(AliStack* const stack, AliESDtrack *const esdTrack)
{
  // Fill DCA comparison information
  if(!esdTrack) return;

  Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
  esdTrack->GetImpactParameters(dca,cov);

  if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
  if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters  
  Int_t clusterITS[200];
  if(esdTrack->GetITSclusters(clusterITS)<fCutsRC->GetMinNClustersITS()) return;  // min. nb. ITS clusters

  Double_t vDCAHisto[4]={dca[0],dca[1],esdTrack->Eta(),esdTrack->Pt()};
  fDCAHisto->Fill(vDCAHisto);

  //
  // Fill rec vs MC information
  //
  if(!stack) return;

}

void AliPerformanceDCA::ProcessConstrained(AliStack* const /*stack*/, AliESDtrack *const /*esdTrack*/)
{
  // Fill DCA comparison information
  
  AliDebug(AliLog::kWarning, "Warning: Not implemented");
}

//_____________________________________________________________________________
Long64_t AliPerformanceDCA::Merge(TCollection* const list) 
{
  // Merge list of objects (needed by PROOF)

  if (!list)
  return 0;

  if (list->IsEmpty())
  return 1;

  TIterator* iter = list->MakeIterator();
  TObject* obj = 0;

  // collection of generated histograms
  Int_t count=0;
  while((obj = iter->Next()) != 0) 
  {
    AliPerformanceDCA* entry = dynamic_cast<AliPerformanceDCA*>(obj);
    if (entry == 0) continue; 

    fDCAHisto->Add(entry->fDCAHisto);
    count++;
  }

return count;
}

//_____________________________________________________________________________
void AliPerformanceDCA::Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, const Bool_t bUseMC)
{
  // Process comparison information
  if(!esdEvent) 
  {
      AliDebug(AliLog::kError, "esdEvent not available");
      return;
  }
  AliHeader* header = 0;
  AliGenEventHeader* genHeader = 0;
  AliStack* stack = 0;
  TArrayF vtxMC(3);

  if(bUseMC)
  {
    if(!mcEvent) {
      AliDebug(AliLog::kError, "mcEvent not available");
      return;
    }
    // get MC event header
    header = mcEvent->Header();
    if (!header) {
      AliDebug(AliLog::kError, "Header not available");
      return;
    }
    // MC particle stack
    stack = mcEvent->Stack();
    if (!stack) {
      AliDebug(AliLog::kError, "Stack not available");
      return;
    }
    // get MC vertex
    genHeader = header->GenEventHeader();
    if (!genHeader) {
      AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
      return;
    }
    genHeader->PrimaryVertex(vtxMC);

  } // end bUseMC

  //  Process events
  for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++) 
  { 
    AliESDtrack *track = esdEvent->GetTrack(iTrack);
    if(!track) continue;

    if(GetAnalysisMode() == 0) ProcessTPC(stack,track);
    else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track);
    else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track);
    else {
      printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
      return;
    }
  }
}

//_____________________________________________________________________________
void AliPerformanceDCA::Analyse()
{
  //
  // Analyse comparison information and store output histograms
  // in the analysis folder "folderDCA" 
  //
  
  TH1::AddDirectory(kFALSE);
  TH1 *h1D=0;
  TH2 *h2D=0;
  //TH3 *h3D=0;
  TObjArray *aFolderObj = new TObjArray;
  char title[256];

  // set pt measurable range 
  fDCAHisto->GetAxis(3)->SetRangeUser(0.10,10.);

  //
  h2D = fDCAHisto->Projection(0,1); // inverse projection convention
  h2D->SetName("dca_r_vs_dca_z");
  h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(1)->GetTitle());
  h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(0)->GetTitle());
  sprintf(title,"%s vs %s",fDCAHisto->GetAxis(0)->GetTitle(),fDCAHisto->GetAxis(1)->GetTitle());
  h2D->SetTitle(title);
  aFolderObj->Add(h2D);

  //
  h2D = fDCAHisto->Projection(0,2);
  h2D->SetName("dca_r_vs_eta");
  h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
  h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(0)->GetTitle());
  sprintf(title,"%s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(0)->GetTitle());
  h2D->SetTitle(title);
  aFolderObj->Add(h2D);

  h1D = MakeStat1D(h2D,0,0);
  h1D->SetName("mean_dca_r_vs_eta");
  h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
  h1D->GetYaxis()->SetTitle("mean_dca_r (cm)");
  sprintf(title," mean_dca_r (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
  h1D->SetTitle(title);
  aFolderObj->Add(h1D);

  h1D = MakeStat1D(h2D,0,1);
  h1D->SetName("rms_dca_r_vs_eta");
  h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
  h1D->GetYaxis()->SetTitle("rms_dca_r (cm)");
  sprintf(title," rms_dca_r (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
  h1D->SetTitle(title);
  aFolderObj->Add(h1D);

  //
  h2D = fDCAHisto->Projection(0,3);
  h2D->SetName("dca_r_vs_pt");
  h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
  h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(0)->GetTitle());
  sprintf(title,"%s vs %s",fDCAHisto->GetAxis(0)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
  h2D->SetTitle(title);
  h2D->SetBit(TH1::kLogX);
  aFolderObj->Add(h2D);

  h1D = MakeStat1D(h2D,0,0);
  h1D->SetName("mean_dca_r_vs_pt");
  h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
  h1D->GetYaxis()->SetTitle("mean_dca_r (cm)");
  sprintf(title,"mean_dca_r (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
  h1D->SetTitle(title);
  h1D->SetBit(TH1::kLogX);
  aFolderObj->Add(h1D);

  h1D = MakeStat1D(h2D,0,1);
  h1D->SetName("rms_dca_r_vs_pt");
  h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
  h1D->GetYaxis()->SetTitle("rms_dca_r (cm)");
  sprintf(title,"rms_dca_r (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
  h1D->SetTitle(title);
  h1D->SetBit(TH1::kLogX);
  aFolderObj->Add(h1D);
   
  // 
  h2D = fDCAHisto->Projection(1,2);
  h2D->SetName("dca_z_vs_eta");
  h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
  h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(1)->GetTitle());
  sprintf(title,"%s vs %s",fDCAHisto->GetAxis(1)->GetTitle(),fDCAHisto->GetAxis(2)->GetTitle());
  h2D->SetTitle(title);
  aFolderObj->Add(h2D);

  h1D = MakeStat1D(h2D,0,0);
  h1D->SetName("mean_dca_z_vs_eta");
  h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
  h1D->GetYaxis()->SetTitle("mean_dca_z (cm)");
  sprintf(title,"mean_dca_z (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
  h1D->SetTitle(title);
  aFolderObj->Add(h1D);

  h1D = MakeStat1D(h2D,0,1);
  h1D->SetName("rms_dca_z_vs_eta");
  h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
  h1D->GetYaxis()->SetTitle("rms_dca_z (cm)");
  sprintf(title,"rms_dca_z (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
  h1D->SetTitle(title);
  aFolderObj->Add(h1D);

  //
  h2D = fDCAHisto->Projection(1,3);
  h2D->SetName("dca_z_vs_pt");
  h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
  h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(1)->GetTitle());
  sprintf(title,"%s vs %s",fDCAHisto->GetAxis(1)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
  h2D->SetTitle(title);
  h2D->SetBit(TH1::kLogX);
  aFolderObj->Add(h2D);

  h1D = MakeStat1D(h2D,0,0);
  h1D->SetName("mean_dca_z_vs_pt");
  h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
  h1D->GetYaxis()->SetTitle("mean_dca_z (cm)");
  sprintf(title,"mean_dca_z (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
  h1D->SetTitle(title);
  h1D->SetBit(TH1::kLogX);
  aFolderObj->Add(h1D);

  h1D = MakeStat1D(h2D,0,1);
  h1D->SetName("rms_dca_z_vs_pt");
  h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
  h1D->GetYaxis()->SetTitle("rms_dca_z (cm)");
  sprintf(title,"rms_dca_z (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
  h1D->SetTitle(title);
  h1D->SetBit(TH1::kLogX);
  aFolderObj->Add(h1D);
   
  /*
  h3D = fDCAHisto->Projection(2,3,0); // normal 3D projection convention
  h3D->SetName("dca_r_vs_eta_vs_pt");

  h2D = MakeStat2D(h3D,0,0,0);
  h2D->SetName("mean_dca_r_vs_eta_vs_pt");
  h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
  h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
  h2D->GetZaxis()->SetTitle("mean_dca_r (cm)");
  sprintf(title,"mean_dca_r (cm) vs %s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
  h2D->SetTitle(title);
  aFolderObj->Add(h2D);

  h2D = MakeStat2D(h3D,0,0,1);
  h2D->SetName("rms_dca_r_vs_eta_vs_pt");
  h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
  h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
  h2D->GetZaxis()->SetTitle("rms_dca_r (cm)");
  sprintf(title,"rms_dca_r (cm) vs %s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
  h2D->SetTitle(title);
  aFolderObj->Add(h2D);

  //
  h3D = fDCAHisto->Projection(2,3,1);
  h3D->SetName("dca_z_vs_eta_vs_pt");

  h2D = MakeStat2D(h3D,0,0,0);
  h2D->SetName("mean_dca_z_vs_eta_vs_pt");
  h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
  h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
  h2D->GetZaxis()->SetTitle("mean_dca_z (cm)");
  sprintf(title,"mean_dca_z (cm) vs %s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
  h2D->SetTitle(title);
  aFolderObj->Add(h2D);

  h2D = MakeStat2D(h3D,0,0,1);
  h2D->SetName("rms_dca_z_vs_eta_vs_pt");
  h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
  h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
  h2D->GetZaxis()->SetTitle("rms_dca_z (cm)");
  sprintf(title,"rms_dca_z (cm) vs %s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
  h2D->SetTitle(title);
  aFolderObj->Add(h2D);
  */


  // export objects to analysis folder
  fAnalysisFolder = ExportToFolder(aFolderObj);

  // delete only TObjArray
  if(aFolderObj) delete aFolderObj;
}

//_____________________________________________________________________________
TH1F* AliPerformanceDCA::MakeStat1D(TH2 *hist, Int_t delta0, Int_t type) 
{
  // Return TH1F histogram 
  // delta - number of bins to integrate
  // with mean (type == 0) or RMS (type==1) 

  char hname[256];
  const char* suffix = "_stat1d";
  sprintf(hname,"%s%s",hist->GetName(),suffix);
  TAxis* xaxis = hist->GetXaxis();
  Int_t  nbinx = xaxis->GetNbins();

  TH1F *hnew = (TH1F*)hist->ProjectionX()->Clone();
  hnew->SetName(hname);

  char name[256];
  for (Int_t ix=0; ix<=nbinx;ix++) {
    sprintf(name,"%s_%d",hist->GetName(),ix);
    TH1 *projection = hist->ProjectionY(name,ix-delta0,ix+delta0);

    Float_t stat= 0., stat_err =0.;
    if (type==0) { stat = projection->GetMean(); stat_err = projection->GetMeanError(); } 
    if (type==1) { stat = projection->GetRMS(); stat_err = projection->GetRMSError(); }
 
    hnew->SetBinContent(ix, stat);
    hnew->SetBinError(ix, stat_err);
  }
  
return hnew;
}

//_____________________________________________________________________________
TH2F* AliPerformanceDCA::MakeStat2D(TH3 *hist, Int_t delta0, Int_t delta1, Int_t type) 
{
  // Return TH1F histogram 
  // delta0 - number of bins to integrate in x
  // delta1 - number of bins to integrate in y
  // with mean (type==0) or RMS (type==1) 

  char hname[256];
  const char* suffix = "_stat2d";
  sprintf(hname,"%s%s",hist->GetName(),suffix);

  TAxis* xaxis = hist->GetXaxis();
  Int_t  nbinx = xaxis->GetNbins(); 

  TH2F *hnew = (TH2F*)hist->Project3D("yx")->Clone();
  hnew->SetName(hname);

  TAxis* yaxis = hist->GetYaxis();
  Int_t  nbiny = yaxis->GetNbins(); 

  char name[256];
  for (Int_t ix=0; ix<=nbinx;ix++) {
    for (Int_t iy=0; iy<=nbiny;iy++) {
      sprintf(name,"%s_%d_%d",hist->GetName(),ix,iy);
      TH1 *projection = hist->ProjectionZ(name,ix-delta0,ix+delta0,iy-delta1,iy+delta1);

      Float_t stat= 0., stat_err =0.;
      if (type==0) { stat = projection->GetMean(); stat_err = projection->GetMeanError(); } 
      if (type==1) { stat = projection->GetRMS(); stat_err = projection->GetRMSError(); }
     
      hnew->SetBinContent(ix,iy,stat);
      hnew->SetBinError(ix,iy,stat_err);
    }
  }
  
return hnew;
}

//_____________________________________________________________________________
TFolder* AliPerformanceDCA::ExportToFolder(TObjArray * array) 
{
  // recreate folder avery time and export objects to new one
  //
  AliPerformanceDCA * comp=this;
  TFolder *folder = comp->GetAnalysisFolder();

  TString name, title;
  TFolder *newFolder = 0;
  Int_t i = 0;
  Int_t size = array->GetSize();

  if(folder) { 
     // get name and title from old folder
     name = folder->GetName();  
     title = folder->GetTitle();  

	 // delete old one
     delete folder;

	 // create new one
     newFolder = CreateFolder(name.Data(),title.Data());
     newFolder->SetOwner();

	 // add objects to folder
     while(i < size) {
	   newFolder->Add(array->At(i));
	   i++;
	 }
  }

return newFolder;
}

//_____________________________________________________________________________
TFolder* AliPerformanceDCA::CreateFolder(TString name,TString title) { 
// create folder for analysed histograms
TFolder *folder = 0;
  folder = new TFolder(name.Data(),title.Data());

  return folder;
}
Commit	Line	Data
777a0ba6	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
	43	#include "AliPerformanceDCA.h"
	44	#include "AliESDEvent.h"
	45	#include "AliESDVertex.h"
	46	#include "AliLog.h"
	47	#include "AliMathBase.h"
	48	#include "AliRecInfoCuts.h"
	49	#include "AliMCInfoCuts.h"
	50	#include "AliStack.h"
	51	#include "AliMCEvent.h"
	52	#include "AliTracker.h"
	53	#include "AliHeader.h"
	54	#include "AliGenEventHeader.h"
	55
	56	using namespace std;
	57
	58	ClassImp(AliPerformanceDCA)
	59
	60	//_____________________________________________________________________________
	61	AliPerformanceDCA::AliPerformanceDCA():
	62	AliPerformanceObject("AliPerformanceDCA"),
	63
	64	// DCA histograms
65	fDCAHisto(0),
66
67	// Cuts
68	fCutsRC(0),
69	fCutsMC(0),
70
71	// histogram folder
72	fAnalysisFolder(0)
73	{
74	// default constructor
75	Init();
76	}
77
78	//_____________________________________________________________________________
79	AliPerformanceDCA::AliPerformanceDCA(Char_t* name="AliPerformanceDCA", Char_t* title="AliPerformanceDCA",Int_t analysisMode=0, Bool_t hptGenerator=kFALSE):
80	AliPerformanceObject(name,title),
81
82	// DCA histograms
83	fDCAHisto(0),
84
85	// Cuts
86	fCutsRC(0),
87	fCutsMC(0),
88
89	// histogram folder
90	fAnalysisFolder(0)
91	{
92	// named constructor
93
94	SetAnalysisMode(analysisMode);
95	SetHptGenerator(hptGenerator);
96	Init();
97	}
98
99	//_____________________________________________________________________________
100	AliPerformanceDCA::~AliPerformanceDCA()
101	{
102	// destructor
103	if(fDCAHisto) delete fDCAHisto; fDCAHisto=0;
104	if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
105	}
106
107	//_____________________________________________________________________________
108	void AliPerformanceDCA::Init()
109	{
110	// DCA histograms
111	Int_t nPtBins = 50;
112	Double_t ptMin = 1.e-2, ptMax = 10.;
113
114	Double_t *binsPt = 0;
115	if (IsHptGenerator()) {
116	nPtBins = 100; ptMax = 100.;
117	binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
118	} else {
119	binsPt = CreateLogAxis(nPtBins,ptMin,ptMax);
120	}
121
122	/*
123	Int_t nPtBins = 31;
124	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.};
125	Double_t ptMin = 0., ptMax = 10.;
126
127	if(IsHptGenerator() == kTRUE) {
128	nPtBins = 100;
129	ptMin = 0.; ptMax = 100.;
130	}
131	*/
132
133	//dca_r, dca_z, eta, pt
134	Int_t binsQA[4] = {100,100,30,nPtBins};
135	Double_t xminQA[4] = {-10.,-10.,-1.5,ptMin};
136	Double_t xmaxQA[4] = {10.,10.,1.5,ptMax};
137
138	fDCAHisto = new THnSparseF("fDCAHisto","dca_r:dca_z:eta:pt",4,binsQA,xminQA,xmaxQA);
139	fDCAHisto->SetBinEdges(3,binsPt);
140
141	fDCAHisto->GetAxis(0)->SetTitle("dca_r (cm)");
142	fDCAHisto->GetAxis(1)->SetTitle("dca_z (cm)");
143	fDCAHisto->GetAxis(2)->SetTitle("#eta");
144	fDCAHisto->GetAxis(3)->SetTitle("p_{T} (GeV/c)");
145	fDCAHisto->Sumw2();
146
147	// init cuts
148	if(!fCutsMC)
149	AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
150	if(!fCutsRC)
151	AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
152
153	// init folder
154	fAnalysisFolder = CreateFolder("folderDCA","Analysis DCA Folder");
155	}
156
157	//_____________________________________________________________________________
158	void AliPerformanceDCA::ProcessTPC(AliStack* const stack, AliESDtrack *const esdTrack)
159	{
160	// Fill DCA comparison information
161	if(!esdTrack) return;
162
163	const AliExternalTrackParam *track = esdTrack->GetTPCInnerParam();
164	if(!track) return;
165
166	Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
167	esdTrack->GetImpactParametersTPC(dca,cov);
168
169	if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
170
171	Double_t vDCAHisto[4]={dca[0],dca[1],track->Eta(),track->Pt()};
172	fDCAHisto->Fill(vDCAHisto);
173
174	//
175	// Fill rec vs MC information
176	//
177	if(!stack) return;
178
179	}
180
181	//_____________________________________________________________________________
182	void AliPerformanceDCA::ProcessTPCITS(AliStack* const stack, AliESDtrack *const esdTrack)
183	{
184	// Fill DCA comparison information
185	if(!esdTrack) return;
186
187	Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
188	esdTrack->GetImpactParameters(dca,cov);
189
190	if ((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
191	if (esdTrack->GetTPCNcls()<fCutsRC->GetMinNClustersTPC()) return; // min. nb. TPC clusters
192	Int_t clusterITS[200];
193	if(esdTrack->GetITSclusters(clusterITS)<fCutsRC->GetMinNClustersITS()) return; // min. nb. ITS clusters
194
195	Double_t vDCAHisto[4]={dca[0],dca[1],esdTrack->Eta(),esdTrack->Pt()};
196	fDCAHisto->Fill(vDCAHisto);
197
198	//
199	// Fill rec vs MC information
200	//
201	if(!stack) return;
202
203	}
204
205	void AliPerformanceDCA::ProcessConstrained(AliStack* const /stack/, AliESDtrack const /esdTrack*/)
206	{
207	// Fill DCA comparison information
208
209	AliDebug(AliLog::kWarning, "Warning: Not implemented");
210	}
211
212	//_____________________________________________________________________________
213	Long64_t AliPerformanceDCA::Merge(TCollection* const list)
214	{
215	// Merge list of objects (needed by PROOF)
216
217	if (!list)
218	return 0;
219
220	if (list->IsEmpty())
221	return 1;
222
223	TIterator* iter = list->MakeIterator();
224	TObject* obj = 0;
225
226	// collection of generated histograms
227	Int_t count=0;
228	while((obj = iter->Next()) != 0)
229	{
230	AliPerformanceDCA* entry = dynamic_cast<AliPerformanceDCA*>(obj);
231	if (entry == 0) continue;
232
233	fDCAHisto->Add(entry->fDCAHisto);
234	count++;
235	}
236
237	return count;
238	}
239
240	//_____________________________________________________________________________
241	void AliPerformanceDCA::Exec(AliMCEvent* const mcEvent, AliESDEvent *const esdEvent, const Bool_t bUseMC)
242	{
243	// Process comparison information
244	if(!esdEvent)
245	{
246	AliDebug(AliLog::kError, "esdEvent not available");
247	return;
248	}
249	AliHeader* header = 0;
250	AliGenEventHeader* genHeader = 0;
251	AliStack* stack = 0;
252	TArrayF vtxMC(3);
253
254	if(bUseMC)
255	{
256	if(!mcEvent) {
257	AliDebug(AliLog::kError, "mcEvent not available");
258	return;
259	}
260	// get MC event header
261	header = mcEvent->Header();
262	if (!header) {
263	AliDebug(AliLog::kError, "Header not available");
264	return;
265	}
266	// MC particle stack
267	stack = mcEvent->Stack();
268	if (!stack) {
269	AliDebug(AliLog::kError, "Stack not available");
270	return;
271	}
272	// get MC vertex
273	genHeader = header->GenEventHeader();
274	if (!genHeader) {
275	AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
276	return;
277	}
278	genHeader->PrimaryVertex(vtxMC);
279
280	} // end bUseMC
281
282	// Process events
283	for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
284	{
285	AliESDtrack *track = esdEvent->GetTrack(iTrack);
286	if(!track) continue;
287
288	if(GetAnalysisMode() == 0) ProcessTPC(stack,track);
289	else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track);
290	else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track);
291	else {
292	printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
293	return;
294	}
295	}
296	}
297
298	//_____________________________________________________________________________
299	void AliPerformanceDCA::Analyse()
300	{
301	//
302	// Analyse comparison information and store output histograms
303	// in the analysis folder "folderDCA"
304	//
305
306	TH1::AddDirectory(kFALSE);
307	TH1 *h1D=0;
308	TH2 *h2D=0;
309	//TH3 *h3D=0;
310	TObjArray *aFolderObj = new TObjArray;
311	char title[256];
312
313	// set pt measurable range
314	fDCAHisto->GetAxis(3)->SetRangeUser(0.10,10.);
315
316	//
317	h2D = fDCAHisto->Projection(0,1); // inverse projection convention
318	h2D->SetName("dca_r_vs_dca_z");
319	h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(1)->GetTitle());
320	h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(0)->GetTitle());
321	sprintf(title,"%s vs %s",fDCAHisto->GetAxis(0)->GetTitle(),fDCAHisto->GetAxis(1)->GetTitle());
322	h2D->SetTitle(title);
323	aFolderObj->Add(h2D);
324
325	//
326	h2D = fDCAHisto->Projection(0,2);
327	h2D->SetName("dca_r_vs_eta");
328	h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
329	h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(0)->GetTitle());
330	sprintf(title,"%s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(0)->GetTitle());
331	h2D->SetTitle(title);
332	aFolderObj->Add(h2D);
333
334	h1D = MakeStat1D(h2D,0,0);
335	h1D->SetName("mean_dca_r_vs_eta");
336	h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
337	h1D->GetYaxis()->SetTitle("mean_dca_r (cm)");
338	sprintf(title," mean_dca_r (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
339	h1D->SetTitle(title);
340	aFolderObj->Add(h1D);
341
342	h1D = MakeStat1D(h2D,0,1);
343	h1D->SetName("rms_dca_r_vs_eta");
344	h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
345	h1D->GetYaxis()->SetTitle("rms_dca_r (cm)");
346	sprintf(title," rms_dca_r (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
347	h1D->SetTitle(title);
348	aFolderObj->Add(h1D);
349
350	//
351	h2D = fDCAHisto->Projection(0,3);
352	h2D->SetName("dca_r_vs_pt");
353	h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
354	h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(0)->GetTitle());
355	sprintf(title,"%s vs %s",fDCAHisto->GetAxis(0)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
356	h2D->SetTitle(title);
357	h2D->SetBit(TH1::kLogX);
358	aFolderObj->Add(h2D);
359
360	h1D = MakeStat1D(h2D,0,0);
361	h1D->SetName("mean_dca_r_vs_pt");
362	h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
363	h1D->GetYaxis()->SetTitle("mean_dca_r (cm)");
364	sprintf(title,"mean_dca_r (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
365	h1D->SetTitle(title);
366	h1D->SetBit(TH1::kLogX);
367	aFolderObj->Add(h1D);
368
369	h1D = MakeStat1D(h2D,0,1);
370	h1D->SetName("rms_dca_r_vs_pt");
371	h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
372	h1D->GetYaxis()->SetTitle("rms_dca_r (cm)");
373	sprintf(title,"rms_dca_r (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
374	h1D->SetTitle(title);
375	h1D->SetBit(TH1::kLogX);
376	aFolderObj->Add(h1D);
377
378	//
379	h2D = fDCAHisto->Projection(1,2);
380	h2D->SetName("dca_z_vs_eta");
381	h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
382	h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(1)->GetTitle());
383	sprintf(title,"%s vs %s",fDCAHisto->GetAxis(1)->GetTitle(),fDCAHisto->GetAxis(2)->GetTitle());
384	h2D->SetTitle(title);
385	aFolderObj->Add(h2D);
386
387	h1D = MakeStat1D(h2D,0,0);
388	h1D->SetName("mean_dca_z_vs_eta");
389	h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
390	h1D->GetYaxis()->SetTitle("mean_dca_z (cm)");
391	sprintf(title,"mean_dca_z (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
392	h1D->SetTitle(title);
393	aFolderObj->Add(h1D);
394
395	h1D = MakeStat1D(h2D,0,1);
396	h1D->SetName("rms_dca_z_vs_eta");
397	h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
398	h1D->GetYaxis()->SetTitle("rms_dca_z (cm)");
399	sprintf(title,"rms_dca_z (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
400	h1D->SetTitle(title);
401	aFolderObj->Add(h1D);
402
403	//
404	h2D = fDCAHisto->Projection(1,3);
405	h2D->SetName("dca_z_vs_pt");
406	h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
407	h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(1)->GetTitle());
408	sprintf(title,"%s vs %s",fDCAHisto->GetAxis(1)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
409	h2D->SetTitle(title);
410	h2D->SetBit(TH1::kLogX);
411	aFolderObj->Add(h2D);
412
413	h1D = MakeStat1D(h2D,0,0);
414	h1D->SetName("mean_dca_z_vs_pt");
415	h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
416	h1D->GetYaxis()->SetTitle("mean_dca_z (cm)");
417	sprintf(title,"mean_dca_z (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
418	h1D->SetTitle(title);
419	h1D->SetBit(TH1::kLogX);
420	aFolderObj->Add(h1D);
421
422	h1D = MakeStat1D(h2D,0,1);
423	h1D->SetName("rms_dca_z_vs_pt");
424	h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
425	h1D->GetYaxis()->SetTitle("rms_dca_z (cm)");
426	sprintf(title,"rms_dca_z (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
427	h1D->SetTitle(title);
428	h1D->SetBit(TH1::kLogX);
429	aFolderObj->Add(h1D);
430
431	/*
432	h3D = fDCAHisto->Projection(2,3,0); // normal 3D projection convention
433	h3D->SetName("dca_r_vs_eta_vs_pt");
434
435	h2D = MakeStat2D(h3D,0,0,0);
436	h2D->SetName("mean_dca_r_vs_eta_vs_pt");
437	h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
438	h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
439	h2D->GetZaxis()->SetTitle("mean_dca_r (cm)");
440	sprintf(title,"mean_dca_r (cm) vs %s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
441	h2D->SetTitle(title);
442	aFolderObj->Add(h2D);
443
444	h2D = MakeStat2D(h3D,0,0,1);
445	h2D->SetName("rms_dca_r_vs_eta_vs_pt");
446	h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
447	h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
448	h2D->GetZaxis()->SetTitle("rms_dca_r (cm)");
449	sprintf(title,"rms_dca_r (cm) vs %s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
450	h2D->SetTitle(title);
451	aFolderObj->Add(h2D);
452
453	//
454	h3D = fDCAHisto->Projection(2,3,1);
455	h3D->SetName("dca_z_vs_eta_vs_pt");
456
457	h2D = MakeStat2D(h3D,0,0,0);
458	h2D->SetName("mean_dca_z_vs_eta_vs_pt");
459	h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
460	h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
461	h2D->GetZaxis()->SetTitle("mean_dca_z (cm)");
462	sprintf(title,"mean_dca_z (cm) vs %s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
463	h2D->SetTitle(title);
464	aFolderObj->Add(h2D);
465
466	h2D = MakeStat2D(h3D,0,0,1);
467	h2D->SetName("rms_dca_z_vs_eta_vs_pt");
468	h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
469	h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
470	h2D->GetZaxis()->SetTitle("rms_dca_z (cm)");
471	sprintf(title,"rms_dca_z (cm) vs %s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
472	h2D->SetTitle(title);
473	aFolderObj->Add(h2D);
474	*/
475
476
477
478
479
480
481
482	// export objects to analysis folder
483	fAnalysisFolder = ExportToFolder(aFolderObj);
484
485	// delete only TObjArray
486	if(aFolderObj) delete aFolderObj;
487	}
488
489	//_____________________________________________________________________________
490	TH1F* AliPerformanceDCA::MakeStat1D(TH2 *hist, Int_t delta0, Int_t type)
491	{
492	// Return TH1F histogram
493	// delta - number of bins to integrate
494	// with mean (type == 0) or RMS (type==1)
495
496	char hname[256];
497	const char* suffix = "_stat1d";
498	sprintf(hname,"%s%s",hist->GetName(),suffix);
499	TAxis* xaxis = hist->GetXaxis();
500	Int_t nbinx = xaxis->GetNbins();
501
502	TH1F hnew = (TH1F)hist->ProjectionX()->Clone();
503	hnew->SetName(hname);
504
505	char name[256];
506	for (Int_t ix=0; ix<=nbinx;ix++) {
507	sprintf(name,"%s_%d",hist->GetName(),ix);
508	TH1 *projection = hist->ProjectionY(name,ix-delta0,ix+delta0);
509
510	Float_t stat= 0., stat_err =0.;
511	if (type==0) { stat = projection->GetMean(); stat_err = projection->GetMeanError(); }
512	if (type==1) { stat = projection->GetRMS(); stat_err = projection->GetRMSError(); }
513
514	hnew->SetBinContent(ix, stat);
515	hnew->SetBinError(ix, stat_err);
516	}
517
518	return hnew;
519	}
520
521	//_____________________________________________________________________________
522	TH2F* AliPerformanceDCA::MakeStat2D(TH3 *hist, Int_t delta0, Int_t delta1, Int_t type)
523	{
524	// Return TH1F histogram
525	// delta0 - number of bins to integrate in x
526	// delta1 - number of bins to integrate in y
527	// with mean (type==0) or RMS (type==1)
528
529	char hname[256];
530	const char* suffix = "_stat2d";
531	sprintf(hname,"%s%s",hist->GetName(),suffix);
532
533	TAxis* xaxis = hist->GetXaxis();
534	Int_t nbinx = xaxis->GetNbins();
535
536	TH2F hnew = (TH2F)hist->Project3D("yx")->Clone();
537	hnew->SetName(hname);
538
539	TAxis* yaxis = hist->GetYaxis();
540	Int_t nbiny = yaxis->GetNbins();
541
542	char name[256];
543	for (Int_t ix=0; ix<=nbinx;ix++) {
544	for (Int_t iy=0; iy<=nbiny;iy++) {
545	sprintf(name,"%s_%d_%d",hist->GetName(),ix,iy);
546	TH1 *projection = hist->ProjectionZ(name,ix-delta0,ix+delta0,iy-delta1,iy+delta1);
547
548	Float_t stat= 0., stat_err =0.;
549	if (type==0) { stat = projection->GetMean(); stat_err = projection->GetMeanError(); }
550	if (type==1) { stat = projection->GetRMS(); stat_err = projection->GetRMSError(); }
551
552	hnew->SetBinContent(ix,iy,stat);
553	hnew->SetBinError(ix,iy,stat_err);
554	}
555	}
556
557	return hnew;
558	}
559
560	//_____________________________________________________________________________
561	TFolder* AliPerformanceDCA::ExportToFolder(TObjArray * array)
562	{
563	// recreate folder avery time and export objects to new one
564	//
565	AliPerformanceDCA * comp=this;
566	TFolder *folder = comp->GetAnalysisFolder();
567
568	TString name, title;
569	TFolder *newFolder = 0;
570	Int_t i = 0;
571	Int_t size = array->GetSize();
572
573	if(folder) {
574	// get name and title from old folder
575	name = folder->GetName();
576	title = folder->GetTitle();
577
578	// delete old one
579	delete folder;
580
581	// create new one
582	newFolder = CreateFolder(name.Data(),title.Data());
583	newFolder->SetOwner();
584
585	// add objects to folder
586	while(i < size) {
587	newFolder->Add(array->At(i));
588	i++;
589	}
590	}
591
592	return newFolder;
593	}
594
595	//_____________________________________________________________________________
596	TFolder* AliPerformanceDCA::CreateFolder(TString name,TString title) {
597	// create folder for analysed histograms
598	TFolder *folder = 0;
599	folder = new TFolder(name.Data(),title.Data());
600
601	return folder;
602	}