[u/mrichter/AliRoot.git] / PWG1 / TPC / AliPerformanceDEdx.cxx

//------------------------------------------------------------------------------
// Implementation of AliPerformanceDEdx 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 of AliPerformanceDEdx.
//  
// 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");
  //AliPerformanceDEdx * compObj = (AliPerformanceDEdx*)f.Get("AliPerformanceDEdx");
  AliPerformanceDEdx * compObj = (AliPerformanceDEdx*)coutput->FindObject("AliPerformanceDEdx");

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

  // the output histograms/graphs will be stored in the folder "folderDEdx" 
  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_DEdx.root"."recreate");
  compObj->Write(); // compObj->GetAnalysisFolder()->Write();
  fout.Close();

*/

#include <TDirectory.h>
#include <TAxis.h>
#include <TCanvas.h>
#include <TH1.h>
#include <TH2.h>
#include <TF1.h>

#include "AliPerformanceDEdx.h" 
#include "AliESDEvent.h"
#include "AliMCEvent.h"
#include "AliESDtrack.h"
#include "AliStack.h"
#include "AliLog.h" 
#include "AliMCInfoCuts.h" 
#include "AliMathBase.h"
#include "AliRecInfoCuts.h" 
#include "AliTreeDraw.h"
#include "AliHeader.h"
#include "AliGenEventHeader.h"

using namespace std;

ClassImp(AliPerformanceDEdx)

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

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

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

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

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

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

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


//_____________________________________________________________________________
AliPerformanceDEdx::~AliPerformanceDEdx()
{
  // destructor
  if(fDeDxHisto)  delete fDeDxHisto; fDeDxHisto=0; 
  if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
}

//_____________________________________________________________________________
void AliPerformanceDEdx::Init()
{
  // Init histograms
  
  // TPC dEdx
  // set pt bins
  Int_t nPBins = 50;
  Double_t pMin = 1.e-2, pMax = 10.;

  Double_t *binsP = 0;
  if (IsHptGenerator())  { 
    nPBins = 100; pMax = 100.;
    binsP = CreateLogAxis(nPBins,pMin,pMax);
  } else {
    binsP = CreateLogAxis(nPBins,pMin,pMax);
  }


  /*
  Int_t nPBins = 31;
    Double_t binsP[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 pMin = 0., pMax = 10.;

    if(IsHptGenerator() == kTRUE) {
      nPBins = 100;
      pMin = 0.; pMax = 100.;
    }
   */

   //dedx:alpha:y:z:snp:tgl:ncls:p
   //dedx:phi:y:z:snp:tgl:ncls:p
   //Int_t binsQA[8]    = {300, 50, 50,  50, 50, 50, 80, nPBins};
   //Double_t xminQA[8] = {0, -4,-20,-250, -1, -2, 0, pMin};
   //Double_t xmaxQA[8] = {300, 4, 20, 250,  1,  2, 160, pMax};
   Int_t binsQA[8]    = {300, 144, 50,  50, 50, 50, 80, nPBins};
   Double_t xminQA[8] = {0, -TMath::Pi(),-20,-250, -1, -2, 0, pMin};
   Double_t xmaxQA[8] = {300, TMath::Pi(), 20, 250,  1,  2, 160, pMax};

   //fDeDxHisto = new THnSparseF("fDeDxHisto","dedx:alpha:y:z:snp:tgl:ncls:momentum",8,binsQA,xminQA,xmaxQA);
   fDeDxHisto = new THnSparseF("fDeDxHisto","dedx:phi:y:z:snp:tgl:ncls:momentum",8,binsQA,xminQA,xmaxQA);
   fDeDxHisto->SetBinEdges(7,binsP);

   fDeDxHisto->GetAxis(0)->SetTitle("dedx (a.u.)");
   fDeDxHisto->GetAxis(1)->SetTitle("#phi (rad)");
   fDeDxHisto->GetAxis(2)->SetTitle("y (cm)");
   fDeDxHisto->GetAxis(3)->SetTitle("z (cm)");
   fDeDxHisto->GetAxis(4)->SetTitle("sin#phi");
   fDeDxHisto->GetAxis(5)->SetTitle("tan#lambda");
   fDeDxHisto->GetAxis(6)->SetTitle("ncls");
   fDeDxHisto->GetAxis(7)->SetTitle("p (GeV/c)");
   fDeDxHisto->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("folderDEdx","Analysis de/dx Folder");
}

//_____________________________________________________________________________
void AliPerformanceDEdx::ProcessTPC(AliStack* const /*stack*/, AliESDtrack *const /*esdTrack*/)
{
  // Fill dE/dx  comparison information
  AliDebug(AliLog::kWarning, "Warning: Not implemented");
}

//_____________________________________________________________________________
void AliPerformanceDEdx::ProcessInnerTPC(AliStack* const stack, AliESDtrack *const esdTrack)
{
  if(!esdTrack) return;

  const AliExternalTrackParam *innerParam =  esdTrack->GetInnerParam();
  if(!innerParam) return;

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

  if((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit

  //
  // select primaries
  //
  Double_t dcaToVertex = -1;
  if( fCutsRC->GetDCAToVertex2D() ) 
  {
      dcaToVertex = TMath::Sqrt(dca[0]*dca[0]/fCutsRC->GetMaxDCAToVertexXY()/fCutsRC->GetMaxDCAToVertexXY()                    + dca[1]*dca[1]/fCutsRC->GetMaxDCAToVertexZ()/fCutsRC->GetMaxDCAToVertexZ()); 
  }
  if(fCutsRC->GetDCAToVertex2D() && dcaToVertex > 1) return;
  if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[0]) > fCutsRC->GetMaxDCAToVertexXY()) return;
  if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[1]) > fCutsRC->GetMaxDCAToVertexZ()) return;

  Float_t dedx = esdTrack->GetTPCsignal();
  Int_t ncls = esdTrack->GetTPCNcls();

  Double_t pt = innerParam->Pt();
  Double_t lam = TMath::ATan2(innerParam->Pz(),innerParam->Pt());
  Double_t p = pt/TMath::Cos(lam);
  //Double_t alpha = innerParam->GetAlpha();
  Double_t phi = TMath::ATan2(innerParam->Py(),innerParam->Px());
  //if(phi<0.) phi += 2.*TMath::Phi();
  Double_t y = innerParam->GetY();
  Double_t z = innerParam->GetZ();
  Double_t snp = innerParam->GetSnp();
  Double_t tgl = innerParam->GetTgl();

  //Double_t vDeDxHisto[8] = {dedx,alpha,y,z,snp,tgl,ncls,p};
  Double_t vDeDxHisto[8] = {dedx,phi,y,z,snp,tgl,ncls,p};
  fDeDxHisto->Fill(vDeDxHisto); 

  if(!stack) return;
}

//_____________________________________________________________________________
void AliPerformanceDEdx::ProcessTPCITS(AliStack* const /*stack*/, AliESDtrack *const /*esdTrack*/)
{
  // Fill dE/dx  comparison information
  
   AliDebug(AliLog::kWarning, "Warning: Not implemented");
}

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

//_____________________________________________________________________________
Long64_t AliPerformanceDEdx::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) 
  {
    AliPerformanceDEdx* entry = dynamic_cast<AliPerformanceDEdx*>(obj);
    if (entry == 0) continue;

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

return count;
}

//_____________________________________________________________________________
void AliPerformanceDEdx::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;
    }
  }

  // trigger
  Bool_t isEventTriggered = esdEvent->IsTriggerClassFired(GetTriggerClass());
  if(!isEventTriggered) return; 

  // get TPC event vertex
  const AliESDVertex *vtxESD = esdEvent->GetPrimaryVertexTPC();
  if(vtxESD && (vtxESD->GetStatus()<=0)) 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 if(GetAnalysisMode() == 3) ProcessInnerTPC(stack,track);
    else {
      printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
      return;
    }
  }
}

//_____________________________________________________________________________
void AliPerformanceDEdx::Analyse()
{
  // Analyze comparison information and store output histograms
  // in the folder "folderDEdx"
  //
  TH1::AddDirectory(kFALSE);
  TH1F *h1D=0;
  TH2F *h2D=0;
  TObjArray *aFolderObj = new TObjArray;

  char name[256];
  char title[256];

  for(Int_t i=1; i<8; i++) { 
    //
    h2D = (TH2F*)fDeDxHisto->Projection(0,i);

    sprintf(name,"h_dedx_%d_vs_%d",0,i);
    h2D->SetName(name);
    sprintf(title,"%s vs %s",fDeDxHisto->GetAxis(0)->GetTitle(),fDeDxHisto->GetAxis(i)->GetTitle());
    h2D->SetTitle(title);
    h2D->GetXaxis()->SetTitle(fDeDxHisto->GetAxis(i)->GetTitle());
    h2D->GetYaxis()->SetTitle(fDeDxHisto->GetAxis(0)->GetTitle());

    if(i==7) h2D->SetBit(TH1::kLogX);
    aFolderObj->Add(h2D);
  }

  // resolution histograms for mips
  //dedx:phi:y:z:snp:tgl:ncls:p
  fDeDxHisto->GetAxis(2)->SetRangeUser(-15.,14.999);
  fDeDxHisto->GetAxis(3)->SetRangeUser(-120.,119.999);
  fDeDxHisto->GetAxis(4)->SetRangeUser(-0.4, 0.399);
  fDeDxHisto->GetAxis(5)->SetRangeUser(-0.9,0.89);
  fDeDxHisto->GetAxis(6)->SetRangeUser(60.,160.);
  fDeDxHisto->GetAxis(7)->SetRangeUser(0.4,0.499);

  h1D=(TH1F*)fDeDxHisto->Projection(0);
  h1D->SetName("dedx_mips");

  h1D->SetTitle("dedx_mips");
  h1D->GetXaxis()->SetTitle(fDeDxHisto->GetAxis(0)->GetTitle());
  aFolderObj->Add(h1D);

  //
  TObjArray *arr[7] = {0};
  TF1 *f1[7] = {0};
  
  for(Int_t i=1; i<8; i++) 
  { 
    arr[i] = new TObjArray;
    f1[i] = new TF1("gaus","gaus");
    //printf("i %d \n",i);

    h2D = (TH2F*)fDeDxHisto->Projection(0,i);

    f1[i]->SetRange(40,60); // should be pion peak
    h2D->FitSlicesY(f1[i],0,-1,10,"QNR",arr[i]); // gaus fit of pion peak

    h1D = (TH1F*)arr[i]->At(1);
    sprintf(name,"mean_dedx_mips_vs_%d",i);
    h1D->SetName(name);
    sprintf(title,"%s vs %s","mean_dedx_mips (a.u.)",fDeDxHisto->GetAxis(i)->GetTitle());
    h1D->SetTitle(title);
    h1D->GetXaxis()->SetTitle(fDeDxHisto->GetAxis(i)->GetTitle());
    h1D->GetYaxis()->SetTitle("mean_dedx_mips (a.u.)");
    //h1D->SetMinimum(40);
    //h1D->SetMaximum(60);

    aFolderObj->Add(h1D);

    h1D = (TH1F*)arr[i]->At(2);
    sprintf(name,"res_dedx_mips_vs_%d",i);
    h1D->SetName(name);
    sprintf(title,"%s vs %s","res_dedx_mips (a.u)",fDeDxHisto->GetAxis(i)->GetTitle());
    h1D->SetTitle(title);
    h1D->GetXaxis()->SetTitle(fDeDxHisto->GetAxis(i)->GetTitle());
    h1D->GetYaxis()->SetTitle("res_dedx_mips (a.u.)");
    //h1D->SetMinimum(0);
    //h1D->SetMaximum(6);

    aFolderObj->Add(h1D);
  }

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

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

//_____________________________________________________________________________
TFolder* AliPerformanceDEdx::ExportToFolder(TObjArray * array) 
{
  // recreate folder avery time and export objects to new one
  //
  AliPerformanceDEdx * 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* AliPerformanceDEdx::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
7cc34f08	1	//------------------------------------------------------------------------------
	2	// Implementation of AliPerformanceDEdx 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 which is
	8	// a data of AliPerformanceDEdx.
	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	//AliPerformanceDEdx * compObj = (AliPerformanceDEdx*)f.Get("AliPerformanceDEdx");
	20	AliPerformanceDEdx * compObj = (AliPerformanceDEdx*)coutput->FindObject("AliPerformanceDEdx");
	21
	22	// Analyse comparison data
	23	compObj->Analyse();
	24
	25	// the output histograms/graphs will be stored in the folder "folderDEdx"
	26	compObj->GetAnalysisFolder()->ls("*");
	27
	28	// user can save whole comparison object (or only folder with anlysed histograms)
	29	// in the seperate output file (e.g.)
	30	TFile fout("Analysed_DEdx.root"."recreate");
	31	compObj->Write(); // compObj->GetAnalysisFolder()->Write();
	32	fout.Close();
	33
	34	*/
	35
	36	#include <TDirectory.h>
	37	#include <TAxis.h>
	38	#include <TCanvas.h>
	39	#include <TH1.h>
	40	#include <TH2.h>
	41	#include <TF1.h>
	42
	43	#include "AliPerformanceDEdx.h"
	44	#include "AliESDEvent.h"
	45	#include "AliMCEvent.h"
	46	#include "AliESDtrack.h"
	47	#include "AliStack.h"
	48	#include "AliLog.h"
	49	#include "AliMCInfoCuts.h"
	50	#include "AliMathBase.h"
	51	#include "AliRecInfoCuts.h"
	52	#include "AliTreeDraw.h"
	53	#include "AliHeader.h"
	54	#include "AliGenEventHeader.h"
	55
	56	using namespace std;
	57
	58	ClassImp(AliPerformanceDEdx)
	59
	60	//_____________________________________________________________________________
	61	AliPerformanceDEdx::AliPerformanceDEdx():
	62	AliPerformanceObject("AliPerformanceDEdx"),
	63
	64	// dEdx
65	fDeDxHisto(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	AliPerformanceDEdx::AliPerformanceDEdx(Char_t* name="AliPerformanceDEdx", Char_t* title="AliPerformanceDEdx",Int_t analysisMode=0, Bool_t hptGenerator=kFALSE):
80	AliPerformanceObject(name,title),
81
82	// dEdx
83	fDeDxHisto(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	//_____________________________________________________________________________
101	AliPerformanceDEdx::~AliPerformanceDEdx()
102	{
103	// destructor
104	if(fDeDxHisto) delete fDeDxHisto; fDeDxHisto=0;
105	if(fAnalysisFolder) delete fAnalysisFolder; fAnalysisFolder=0;
106	}
107
108	//_____________________________________________________________________________
109	void AliPerformanceDEdx::Init()
110	{
111	// Init histograms
112
113	// TPC dEdx
114	// set pt bins
115	Int_t nPBins = 50;
116	Double_t pMin = 1.e-2, pMax = 10.;
117
118	Double_t *binsP = 0;
119	if (IsHptGenerator()) {
120	nPBins = 100; pMax = 100.;
121	binsP = CreateLogAxis(nPBins,pMin,pMax);
122	} else {
123	binsP = CreateLogAxis(nPBins,pMin,pMax);
124	}
125
126
127	/*
128	Int_t nPBins = 31;
129	Double_t binsP[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.};
130	Double_t pMin = 0., pMax = 10.;
131
132	if(IsHptGenerator() == kTRUE) {
133	nPBins = 100;
134	pMin = 0.; pMax = 100.;
135	}
136	*/
137
138	//dedx:alpha:y:z:snp:tgl:ncls:p
139	//dedx:phi:y:z:snp:tgl:ncls:p
140	//Int_t binsQA[8] = {300, 50, 50, 50, 50, 50, 80, nPBins};
141	//Double_t xminQA[8] = {0, -4,-20,-250, -1, -2, 0, pMin};
142	//Double_t xmaxQA[8] = {300, 4, 20, 250, 1, 2, 160, pMax};
e6a60a90	143	Int_t binsQA[8] = {300, 144, 50, 50, 50, 50, 80, nPBins};
7cc34f08	144	Double_t xminQA[8] = {0, -TMath::Pi(),-20,-250, -1, -2, 0, pMin};
	145	Double_t xmaxQA[8] = {300, TMath::Pi(), 20, 250, 1, 2, 160, pMax};
	146
	147	//fDeDxHisto = new THnSparseF("fDeDxHisto","dedx:alpha:y:z:snp:tgl:ncls:momentum",8,binsQA,xminQA,xmaxQA);
	148	fDeDxHisto = new THnSparseF("fDeDxHisto","dedx:phi:y:z:snp:tgl:ncls:momentum",8,binsQA,xminQA,xmaxQA);
	149	fDeDxHisto->SetBinEdges(7,binsP);
	150
	151	fDeDxHisto->GetAxis(0)->SetTitle("dedx (a.u.)");
	152	fDeDxHisto->GetAxis(1)->SetTitle("#phi (rad)");
	153	fDeDxHisto->GetAxis(2)->SetTitle("y (cm)");
	154	fDeDxHisto->GetAxis(3)->SetTitle("z (cm)");
	155	fDeDxHisto->GetAxis(4)->SetTitle("sin#phi");
	156	fDeDxHisto->GetAxis(5)->SetTitle("tan#lambda");
	157	fDeDxHisto->GetAxis(6)->SetTitle("ncls");
	158	fDeDxHisto->GetAxis(7)->SetTitle("p (GeV/c)");
	159	fDeDxHisto->Sumw2();
	160
	161	// Init cuts
	162	if(!fCutsMC)
	163	AliDebug(AliLog::kError, "ERROR: Cannot find AliMCInfoCuts object");
	164	if(!fCutsRC)
	165	AliDebug(AliLog::kError, "ERROR: Cannot find AliRecInfoCuts object");
	166
	167	// init folder
	168	fAnalysisFolder = CreateFolder("folderDEdx","Analysis de/dx Folder");
	169	}
	170
	171	//_____________________________________________________________________________
	172	void AliPerformanceDEdx::ProcessTPC(AliStack* const /stack/, AliESDtrack const /esdTrack*/)
	173	{
	174	// Fill dE/dx comparison information
	175	AliDebug(AliLog::kWarning, "Warning: Not implemented");
	176	}
	177
	178	//_____________________________________________________________________________
	179	void AliPerformanceDEdx::ProcessInnerTPC(AliStack* const stack, AliESDtrack *const esdTrack)
	180	{
	181	if(!esdTrack) return;
	182
	183	const AliExternalTrackParam *innerParam = esdTrack->GetInnerParam();
	184	if(!innerParam) return;
	185
	186	Float_t dca[2], cov[3]; // dca_xy, dca_z, sigma_xy, sigma_xy_z, sigma_z
	187	esdTrack->GetImpactParametersTPC(dca,cov);
	188
	189	if((esdTrack->GetStatus()&AliESDtrack::kTPCrefit)==0) return; // TPC refit
a13f4653	190
	191	//
	192	// select primaries
	193	//
	194	Double_t dcaToVertex = -1;
	195	if( fCutsRC->GetDCAToVertex2D() )
	196	{
	197	dcaToVertex = TMath::Sqrt(dca[0]dca[0]/fCutsRC->GetMaxDCAToVertexXY()/fCutsRC->GetMaxDCAToVertexXY() + dca[1]dca[1]/fCutsRC->GetMaxDCAToVertexZ()/fCutsRC->GetMaxDCAToVertexZ());
7cc34f08	198	}
a13f4653	199	if(fCutsRC->GetDCAToVertex2D() && dcaToVertex > 1) return;
	200	if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[0]) > fCutsRC->GetMaxDCAToVertexXY()) return;
	201	if(!fCutsRC->GetDCAToVertex2D() && TMath::Abs(dca[1]) > fCutsRC->GetMaxDCAToVertexZ()) return;
	202
	203	Float_t dedx = esdTrack->GetTPCsignal();
	204	Int_t ncls = esdTrack->GetTPCNcls();
	205
	206	Double_t pt = innerParam->Pt();
	207	Double_t lam = TMath::ATan2(innerParam->Pz(),innerParam->Pt());
	208	Double_t p = pt/TMath::Cos(lam);
	209	//Double_t alpha = innerParam->GetAlpha();
	210	Double_t phi = TMath::ATan2(innerParam->Py(),innerParam->Px());
	211	//if(phi<0.) phi += 2.*TMath::Phi();
	212	Double_t y = innerParam->GetY();
	213	Double_t z = innerParam->GetZ();
	214	Double_t snp = innerParam->GetSnp();
	215	Double_t tgl = innerParam->GetTgl();
	216
	217	//Double_t vDeDxHisto[8] = {dedx,alpha,y,z,snp,tgl,ncls,p};
	218	Double_t vDeDxHisto[8] = {dedx,phi,y,z,snp,tgl,ncls,p};
	219	fDeDxHisto->Fill(vDeDxHisto);
7cc34f08	220
	221	if(!stack) return;
	222	}
	223
	224	//_____________________________________________________________________________
	225	void AliPerformanceDEdx::ProcessTPCITS(AliStack* const /stack/, AliESDtrack const /esdTrack*/)
	226	{
	227	// Fill dE/dx comparison information
	228
	229	AliDebug(AliLog::kWarning, "Warning: Not implemented");
	230	}
	231
	232	//_____________________________________________________________________________
	233	void AliPerformanceDEdx::ProcessConstrained(AliStack* const /stack/, AliESDtrack const /esdTrack*/)
	234	{
	235	// Fill dE/dx comparison information
	236
	237	AliDebug(AliLog::kWarning, "Warning: Not implemented");
	238	}
	239
	240	//_____________________________________________________________________________
	241	Long64_t AliPerformanceDEdx::Merge(TCollection* const list)
	242	{
	243	// Merge list of objects (needed by PROOF)
	244
	245	if (!list)
	246	return 0;
	247
	248	if (list->IsEmpty())
	249	return 1;
	250
	251	TIterator* iter = list->MakeIterator();
	252	TObject* obj = 0;
	253
	254	// collection of generated histograms
	255	Int_t count=0;
	256	while((obj = iter->Next()) != 0)
	257	{
	258	AliPerformanceDEdx* entry = dynamic_cast<AliPerformanceDEdx*>(obj);
	259	if (entry == 0) continue;
	260
	261	fDeDxHisto->Add(entry->fDeDxHisto);
	262	count++;
	263	}
	264
	265	return count;
	266	}
	267
	268	//_____________________________________________________________________________
	269	void AliPerformanceDEdx::Exec(AliMCEvent* const mcEvent, AliESDEvent const esdEvent, AliESDfriend const esdFriend, const Bool_t bUseMC, const Bool_t bUseESDfriend)
	270	{
	271	// Process comparison information
	272	//
	273	if(!esdEvent)
	274	{
	275	AliDebug(AliLog::kError, "esdEvent not available");
	276	return;
	277	}
	278	AliHeader* header = 0;
	279	AliGenEventHeader* genHeader = 0;
	280	AliStack* stack = 0;
	281	TArrayF vtxMC(3);
	282
	283	if(bUseMC)
284	{
285	if(!mcEvent) {
286	AliDebug(AliLog::kError, "mcEvent not available");
287	return;
288	}
289
290	// get MC event header
291	header = mcEvent->Header();
292	if (!header) {
293	AliDebug(AliLog::kError, "Header not available");
294	return;
295	}
296	// MC particle stack
297	stack = mcEvent->Stack();
298	if (!stack) {
299	AliDebug(AliLog::kError, "Stack not available");
300	return;
301	}
302
303	// get MC vertex
304	genHeader = header->GenEventHeader();
305	if (!genHeader) {
306	AliDebug(AliLog::kError, "Could not retrieve genHeader from Header");
307	return;
308	}
309	genHeader->PrimaryVertex(vtxMC);
310
311	} // end bUseMC
312
313	// use ESD friends
314	if(bUseESDfriend) {
315	if(!esdFriend) {
316	AliDebug(AliLog::kError, "esdFriend not available");
317	return;
318	}
319	}
320
e6a60a90	321	// trigger
	322	Bool_t isEventTriggered = esdEvent->IsTriggerClassFired(GetTriggerClass());
	323	if(!isEventTriggered) return;
7cc34f08	324
e6a60a90	325	// get TPC event vertex
	326	const AliESDVertex *vtxESD = esdEvent->GetPrimaryVertexTPC();
	327	if(vtxESD && (vtxESD->GetStatus()<=0)) return;
7cc34f08	328
	329	// Process events
	330	for (Int_t iTrack = 0; iTrack < esdEvent->GetNumberOfTracks(); iTrack++)
	331	{
	332	AliESDtrack *track = esdEvent->GetTrack(iTrack);
	333	if(!track) continue;
	334
	335	if(GetAnalysisMode() == 0) ProcessTPC(stack,track);
	336	else if(GetAnalysisMode() == 1) ProcessTPCITS(stack,track);
	337	else if(GetAnalysisMode() == 2) ProcessConstrained(stack,track);
	338	else if(GetAnalysisMode() == 3) ProcessInnerTPC(stack,track);
	339	else {
	340	printf("ERROR: AnalysisMode %d \n",fAnalysisMode);
	341	return;
	342	}
	343	}
	344	}
	345
	346	//_____________________________________________________________________________
	347	void AliPerformanceDEdx::Analyse()
	348	{
	349	// Analyze comparison information and store output histograms
	350	// in the folder "folderDEdx"
	351	//
	352	TH1::AddDirectory(kFALSE);
	353	TH1F *h1D=0;
	354	TH2F *h2D=0;
	355	TObjArray *aFolderObj = new TObjArray;
	356
	357	char name[256];
	358	char title[256];
	359
	360	for(Int_t i=1; i<8; i++) {
	361	//
	362	h2D = (TH2F*)fDeDxHisto->Projection(0,i);
	363
	364	sprintf(name,"h_dedx_%d_vs_%d",0,i);
	365	h2D->SetName(name);
	366	sprintf(title,"%s vs %s",fDeDxHisto->GetAxis(0)->GetTitle(),fDeDxHisto->GetAxis(i)->GetTitle());
	367	h2D->SetTitle(title);
	368	h2D->GetXaxis()->SetTitle(fDeDxHisto->GetAxis(i)->GetTitle());
	369	h2D->GetYaxis()->SetTitle(fDeDxHisto->GetAxis(0)->GetTitle());
	370
	371	if(i==7) h2D->SetBit(TH1::kLogX);
	372	aFolderObj->Add(h2D);
	373	}
	374
	375	// resolution histograms for mips
a13f4653	376	//dedx:phi:y:z:snp:tgl:ncls:p
7cc34f08	377	fDeDxHisto->GetAxis(2)->SetRangeUser(-15.,14.999);
	378	fDeDxHisto->GetAxis(3)->SetRangeUser(-120.,119.999);
	379	fDeDxHisto->GetAxis(4)->SetRangeUser(-0.4, 0.399);
a13f4653	380	fDeDxHisto->GetAxis(5)->SetRangeUser(-0.9,0.89);
a13f4653	381	fDeDxHisto->GetAxis(6)->SetRangeUser(60.,160.);
7cc34f08	382	fDeDxHisto->GetAxis(7)->SetRangeUser(0.4,0.499);
	383
	384	h1D=(TH1F*)fDeDxHisto->Projection(0);
	385	h1D->SetName("dedx_mips");
	386
	387	h1D->SetTitle("dedx_mips");
	388	h1D->GetXaxis()->SetTitle(fDeDxHisto->GetAxis(0)->GetTitle());
	389	aFolderObj->Add(h1D);
	390
	391	//
	392	TObjArray *arr[7] = {0};
	393	TF1 *f1[7] = {0};
	394
	395	for(Int_t i=1; i<8; i++)
	396	{
	397	arr[i] = new TObjArray;
	398	f1[i] = new TF1("gaus","gaus");
	399	//printf("i %d \n",i);
	400
	401	h2D = (TH2F*)fDeDxHisto->Projection(0,i);
	402
	403	f1[i]->SetRange(40,60); // should be pion peak
	404	h2D->FitSlicesY(f1[i],0,-1,10,"QNR",arr[i]); // gaus fit of pion peak
	405
	406	h1D = (TH1F*)arr[i]->At(1);
	407	sprintf(name,"mean_dedx_mips_vs_%d",i);
	408	h1D->SetName(name);
	409	sprintf(title,"%s vs %s","mean_dedx_mips (a.u.)",fDeDxHisto->GetAxis(i)->GetTitle());
	410	h1D->SetTitle(title);
	411	h1D->GetXaxis()->SetTitle(fDeDxHisto->GetAxis(i)->GetTitle());
	412	h1D->GetYaxis()->SetTitle("mean_dedx_mips (a.u.)");
	413	//h1D->SetMinimum(40);
	414	//h1D->SetMaximum(60);
	415
	416	aFolderObj->Add(h1D);
	417
	418	h1D = (TH1F*)arr[i]->At(2);
	419	sprintf(name,"res_dedx_mips_vs_%d",i);
	420	h1D->SetName(name);
	421	sprintf(title,"%s vs %s","res_dedx_mips (a.u)",fDeDxHisto->GetAxis(i)->GetTitle());
	422	h1D->SetTitle(title);
	423	h1D->GetXaxis()->SetTitle(fDeDxHisto->GetAxis(i)->GetTitle());
	424	h1D->GetYaxis()->SetTitle("res_dedx_mips (a.u.)");
	425	//h1D->SetMinimum(0);
	426	//h1D->SetMaximum(6);
	427
	428	aFolderObj->Add(h1D);
	429	}
	430
	431	// export objects to analysis folder
	432	fAnalysisFolder = ExportToFolder(aFolderObj);
	433
a13f4653	434	// delete only TObjArray
7cc34f08	435	if(aFolderObj) delete aFolderObj;
	436	}
	437
	438	//_____________________________________________________________________________
	439	TFolder* AliPerformanceDEdx::ExportToFolder(TObjArray * array)
	440	{
	441	// recreate folder avery time and export objects to new one
	442	//
	443	AliPerformanceDEdx * comp=this;
	444	TFolder *folder = comp->GetAnalysisFolder();
	445
	446	TString name, title;
	447	TFolder *newFolder = 0;
	448	Int_t i = 0;
	449	Int_t size = array->GetSize();
	450
	451	if(folder) {
	452	// get name and title from old folder
	453	name = folder->GetName();
	454	title = folder->GetTitle();
	455
	456	// delete old one
	457	delete folder;
	458
	459	// create new one
	460	newFolder = CreateFolder(name.Data(),title.Data());
	461	newFolder->SetOwner();
	462
	463	// add objects to folder
	464	while(i < size) {
	465	newFolder->Add(array->At(i));
	466	i++;
	467	}
	468	}
	469
	470	return newFolder;
	471	}
	472
	473
	474	//_____________________________________________________________________________
	475	TFolder* AliPerformanceDEdx::CreateFolder(TString name,TString title) {
	476	// create folder for analysed histograms
	477	TFolder *folder = 0;
	478	folder = new TFolder(name.Data(),title.Data());
	479
	480	return folder;
	481	}
	482