[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, AliESDfriend *const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend)
{
  // 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

  // use ESD friends
  if(bUseESDfriend) {
    if(!esdFriend) {
      AliDebug(AliLog::kError, "esdFriend not available");
      return;
    }
  }

  //  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	//_____________________________________________________________________________
4f96d707	241	void AliPerformanceDCA::Exec(AliMCEvent* const mcEvent, AliESDEvent const esdEvent, AliESDfriend const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend)
777a0ba6	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
4f96d707	282	// use ESD friends
	283	if(bUseESDfriend) {
	284	if(!esdFriend) {
	285	AliDebug(AliLog::kError, "esdFriend not available");
	286	return;
	287	}
	288	}
	289
777a0ba6	290	// Process events
	291	for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
	292	{
	293	AliESDtrack *track = esdEvent->GetTrack(iTrack);
	294	if(!track) continue;
	295
	296	if(GetAnalysisMode() == 0) ProcessTPC(stack,track);
	297	else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track);
	298	else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track);
	299	else {
	300	printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
	301	return;
	302	}
	303	}
	304	}
	305
	306	//_____________________________________________________________________________
	307	void AliPerformanceDCA::Analyse()
	308	{
	309	//
	310	// Analyse comparison information and store output histograms
	311	// in the analysis folder "folderDCA"
	312	//
	313
	314	TH1::AddDirectory(kFALSE);
	315	TH1 *h1D=0;
	316	TH2 *h2D=0;
	317	//TH3 *h3D=0;
	318	TObjArray *aFolderObj = new TObjArray;
	319	char title[256];
	320
	321	// set pt measurable range
	322	fDCAHisto->GetAxis(3)->SetRangeUser(0.10,10.);
	323
	324	//
	325	h2D = fDCAHisto->Projection(0,1); // inverse projection convention
	326	h2D->SetName("dca_r_vs_dca_z");
	327	h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(1)->GetTitle());
	328	h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(0)->GetTitle());
	329	sprintf(title,"%s vs %s",fDCAHisto->GetAxis(0)->GetTitle(),fDCAHisto->GetAxis(1)->GetTitle());
	330	h2D->SetTitle(title);
	331	aFolderObj->Add(h2D);
	332
	333	//
	334	h2D = fDCAHisto->Projection(0,2);
	335	h2D->SetName("dca_r_vs_eta");
	336	h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
	337	h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(0)->GetTitle());
	338	sprintf(title,"%s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(0)->GetTitle());
	339	h2D->SetTitle(title);
	340	aFolderObj->Add(h2D);
	341
	342	h1D = MakeStat1D(h2D,0,0);
	343	h1D->SetName("mean_dca_r_vs_eta");
	344	h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
	345	h1D->GetYaxis()->SetTitle("mean_dca_r (cm)");
	346	sprintf(title," mean_dca_r (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
	347	h1D->SetTitle(title);
	348	aFolderObj->Add(h1D);
	349
	350	h1D = MakeStat1D(h2D,0,1);
	351	h1D->SetName("rms_dca_r_vs_eta");
	352	h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
	353	h1D->GetYaxis()->SetTitle("rms_dca_r (cm)");
354	sprintf(title," rms_dca_r (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
355	h1D->SetTitle(title);
356	aFolderObj->Add(h1D);
357
358	//
359	h2D = fDCAHisto->Projection(0,3);
360	h2D->SetName("dca_r_vs_pt");
361	h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
362	h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(0)->GetTitle());
363	sprintf(title,"%s vs %s",fDCAHisto->GetAxis(0)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
364	h2D->SetTitle(title);
365	h2D->SetBit(TH1::kLogX);
366	aFolderObj->Add(h2D);
367
368	h1D = MakeStat1D(h2D,0,0);
369	h1D->SetName("mean_dca_r_vs_pt");
370	h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
371	h1D->GetYaxis()->SetTitle("mean_dca_r (cm)");
372	sprintf(title,"mean_dca_r (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
373	h1D->SetTitle(title);
374	h1D->SetBit(TH1::kLogX);
375	aFolderObj->Add(h1D);
376
377	h1D = MakeStat1D(h2D,0,1);
378	h1D->SetName("rms_dca_r_vs_pt");
379	h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
380	h1D->GetYaxis()->SetTitle("rms_dca_r (cm)");
381	sprintf(title,"rms_dca_r (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
382	h1D->SetTitle(title);
383	h1D->SetBit(TH1::kLogX);
384	aFolderObj->Add(h1D);
385
386	//
387	h2D = fDCAHisto->Projection(1,2);
388	h2D->SetName("dca_z_vs_eta");
389	h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
390	h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(1)->GetTitle());
391	sprintf(title,"%s vs %s",fDCAHisto->GetAxis(1)->GetTitle(),fDCAHisto->GetAxis(2)->GetTitle());
392	h2D->SetTitle(title);
393	aFolderObj->Add(h2D);
394
395	h1D = MakeStat1D(h2D,0,0);
396	h1D->SetName("mean_dca_z_vs_eta");
397	h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
398	h1D->GetYaxis()->SetTitle("mean_dca_z (cm)");
399	sprintf(title,"mean_dca_z (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
400	h1D->SetTitle(title);
401	aFolderObj->Add(h1D);
402
403	h1D = MakeStat1D(h2D,0,1);
404	h1D->SetName("rms_dca_z_vs_eta");
405	h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
406	h1D->GetYaxis()->SetTitle("rms_dca_z (cm)");
407	sprintf(title,"rms_dca_z (cm) vs %s",fDCAHisto->GetAxis(2)->GetTitle());
408	h1D->SetTitle(title);
409	aFolderObj->Add(h1D);
410
411	//
412	h2D = fDCAHisto->Projection(1,3);
413	h2D->SetName("dca_z_vs_pt");
414	h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
415	h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(1)->GetTitle());
416	sprintf(title,"%s vs %s",fDCAHisto->GetAxis(1)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
417	h2D->SetTitle(title);
418	h2D->SetBit(TH1::kLogX);
419	aFolderObj->Add(h2D);
420
421	h1D = MakeStat1D(h2D,0,0);
422	h1D->SetName("mean_dca_z_vs_pt");
423	h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
424	h1D->GetYaxis()->SetTitle("mean_dca_z (cm)");
425	sprintf(title,"mean_dca_z (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
426	h1D->SetTitle(title);
427	h1D->SetBit(TH1::kLogX);
428	aFolderObj->Add(h1D);
429
430	h1D = MakeStat1D(h2D,0,1);
431	h1D->SetName("rms_dca_z_vs_pt");
432	h1D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
433	h1D->GetYaxis()->SetTitle("rms_dca_z (cm)");
434	sprintf(title,"rms_dca_z (cm) vs %s",fDCAHisto->GetAxis(3)->GetTitle());
435	h1D->SetTitle(title);
436	h1D->SetBit(TH1::kLogX);
437	aFolderObj->Add(h1D);
438
439	/*
440	h3D = fDCAHisto->Projection(2,3,0); // normal 3D projection convention
441	h3D->SetName("dca_r_vs_eta_vs_pt");
442
443	h2D = MakeStat2D(h3D,0,0,0);
444	h2D->SetName("mean_dca_r_vs_eta_vs_pt");
445	h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
446	h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
447	h2D->GetZaxis()->SetTitle("mean_dca_r (cm)");
448	sprintf(title,"mean_dca_r (cm) vs %s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
449	h2D->SetTitle(title);
450	aFolderObj->Add(h2D);
451
452	h2D = MakeStat2D(h3D,0,0,1);
453	h2D->SetName("rms_dca_r_vs_eta_vs_pt");
454	h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
455	h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
456	h2D->GetZaxis()->SetTitle("rms_dca_r (cm)");
457	sprintf(title,"rms_dca_r (cm) vs %s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
458	h2D->SetTitle(title);
459	aFolderObj->Add(h2D);
460
461	//
462	h3D = fDCAHisto->Projection(2,3,1);
463	h3D->SetName("dca_z_vs_eta_vs_pt");
464
465	h2D = MakeStat2D(h3D,0,0,0);
466	h2D->SetName("mean_dca_z_vs_eta_vs_pt");
467	h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
468	h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
469	h2D->GetZaxis()->SetTitle("mean_dca_z (cm)");
470	sprintf(title,"mean_dca_z (cm) vs %s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
471	h2D->SetTitle(title);
472	aFolderObj->Add(h2D);
473
474	h2D = MakeStat2D(h3D,0,0,1);
475	h2D->SetName("rms_dca_z_vs_eta_vs_pt");
476	h2D->GetXaxis()->SetTitle(fDCAHisto->GetAxis(2)->GetTitle());
477	h2D->GetYaxis()->SetTitle(fDCAHisto->GetAxis(3)->GetTitle());
478	h2D->GetZaxis()->SetTitle("rms_dca_z (cm)");
479	sprintf(title,"rms_dca_z (cm) vs %s vs %s",fDCAHisto->GetAxis(2)->GetTitle(),fDCAHisto->GetAxis(3)->GetTitle());
480	h2D->SetTitle(title);
481	aFolderObj->Add(h2D);
482	*/
483
484
485
486
487
488
489
490	// export objects to analysis folder
491	fAnalysisFolder = ExportToFolder(aFolderObj);
492
493	// delete only TObjArray
494	if(aFolderObj) delete aFolderObj;
495	}
496
497	//_____________________________________________________________________________
498	TH1F* AliPerformanceDCA::MakeStat1D(TH2 *hist, Int_t delta0, Int_t type)
499	{
500	// Return TH1F histogram
501	// delta - number of bins to integrate
502	// with mean (type == 0) or RMS (type==1)
503
504	char hname[256];
505	const char* suffix = "_stat1d";
506	sprintf(hname,"%s%s",hist->GetName(),suffix);
507	TAxis* xaxis = hist->GetXaxis();
508	Int_t nbinx = xaxis->GetNbins();
509
510	TH1F hnew = (TH1F)hist->ProjectionX()->Clone();
511	hnew->SetName(hname);
512
513	char name[256];
514	for (Int_t ix=0; ix<=nbinx;ix++) {
515	sprintf(name,"%s_%d",hist->GetName(),ix);
516	TH1 *projection = hist->ProjectionY(name,ix-delta0,ix+delta0);
517
518	Float_t stat= 0., stat_err =0.;
519	if (type==0) { stat = projection->GetMean(); stat_err = projection->GetMeanError(); }
520	if (type==1) { stat = projection->GetRMS(); stat_err = projection->GetRMSError(); }
521
522	hnew->SetBinContent(ix, stat);
523	hnew->SetBinError(ix, stat_err);
524	}
525
526	return hnew;
527	}
528
529	//_____________________________________________________________________________
530	TH2F* AliPerformanceDCA::MakeStat2D(TH3 *hist, Int_t delta0, Int_t delta1, Int_t type)
531	{
532	// Return TH1F histogram
533	// delta0 - number of bins to integrate in x
534	// delta1 - number of bins to integrate in y
535	// with mean (type==0) or RMS (type==1)
536
537	char hname[256];
538	const char* suffix = "_stat2d";
539	sprintf(hname,"%s%s",hist->GetName(),suffix);
540
541	TAxis* xaxis = hist->GetXaxis();
542	Int_t nbinx = xaxis->GetNbins();
543
544	TH2F hnew = (TH2F)hist->Project3D("yx")->Clone();
545	hnew->SetName(hname);
546
547	TAxis* yaxis = hist->GetYaxis();
548	Int_t nbiny = yaxis->GetNbins();
549
550	char name[256];
551	for (Int_t ix=0; ix<=nbinx;ix++) {
552	for (Int_t iy=0; iy<=nbiny;iy++) {
553	sprintf(name,"%s_%d_%d",hist->GetName(),ix,iy);
554	TH1 *projection = hist->ProjectionZ(name,ix-delta0,ix+delta0,iy-delta1,iy+delta1);
555
556	Float_t stat= 0., stat_err =0.;
557	if (type==0) { stat = projection->GetMean(); stat_err = projection->GetMeanError(); }
558	if (type==1) { stat = projection->GetRMS(); stat_err = projection->GetRMSError(); }
559
560	hnew->SetBinContent(ix,iy,stat);
561	hnew->SetBinError(ix,iy,stat_err);
562	}
563	}
564
565	return hnew;
566	}
567
568	//_____________________________________________________________________________
569	TFolder* AliPerformanceDCA::ExportToFolder(TObjArray * array)
570	{
571	// recreate folder avery time and export objects to new one
572	//
573	AliPerformanceDCA * comp=this;
574	TFolder *folder = comp->GetAnalysisFolder();
575
576	TString name, title;
577	TFolder *newFolder = 0;
578	Int_t i = 0;
579	Int_t size = array->GetSize();
580
581	if(folder) {
582	// get name and title from old folder
583	name = folder->GetName();
584	title = folder->GetTitle();
585
586	// delete old one
587	delete folder;
588
589	// create new one
590	newFolder = CreateFolder(name.Data(),title.Data());
591	newFolder->SetOwner();
592
593	// add objects to folder
594	while(i < size) {
595	newFolder->Add(array->At(i));
596	i++;
597	}
598	}
599
600	return newFolder;
601	}
602
603	//_____________________________________________________________________________
604	TFolder* AliPerformanceDCA::CreateFolder(TString name,TString title) {
605	// create folder for analysed histograms
606	TFolder *folder = 0;
607	folder = new TFolder(name.Data(),title.Data());
608
609	return folder;
610	}