[u/mrichter/AliRoot.git] / FMD / analysis / AliFMDAnalysisTaskCollector.cxx

 
#include <TROOT.h>
#include <TSystem.h>
#include <TInterpreter.h>
#include <TChain.h>
#include <TFile.h>
#include <TList.h>
#include <TMath.h>
#include <iostream>

#include "AliFMDAnalysisTaskCollector.h"
#include "AliAnalysisManager.h"
#include "AliESDFMD.h"
#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliAODHandler.h"
#include "AliMCEventHandler.h"
#include "AliStack.h"
#include "AliESDVertex.h"
#include "AliFMDAnaParameters.h"
//#include "AliFMDGeometry.h"


ClassImp(AliFMDAnalysisTaskCollector)

//____________________________________________________________________
Double_t  AliFMDAnalysisTaskCollector::TripleLandau(Double_t *x, Double_t *par) {
  
  Double_t energy        = x[0];
  Double_t constant = par[0];
  Double_t mpv      = par[1];
  Double_t sigma    = par[2];
  Double_t alpha    = par[3];
  Double_t beta     = par[4];
 
  Double_t f = constant*(TMath::Landau(energy,mpv,sigma,kTRUE)+
			 alpha*TMath::Landau(energy,2*mpv+2*sigma*TMath::Log(2),2*sigma,kTRUE)+
			 beta*TMath::Landau(energy,3*mpv+3*sigma*TMath::Log(3),3*sigma,kTRUE) );
  
  return f;
}
//____________________________________________________________________

AliFMDAnalysisTaskCollector::AliFMDAnalysisTaskCollector()
: fDebug(0),
  fOutputList(0),
  fArray(0),
  fZvtxDist(0)
{
  // Default constructor
  
 
}
//____________________________________________________________________
AliFMDAnalysisTaskCollector::AliFMDAnalysisTaskCollector(const char* name):
    AliAnalysisTaskSE(name),
    fDebug(0),
    fOutputList(0),
    fArray(0),
    fZvtxDist(0)
{
  // Default constructor
  
  DefineOutput(1, TList::Class());
}
//____________________________________________________________________
void AliFMDAnalysisTaskCollector::UserCreateOutputObjects()
{
  // Create the output container
  printf("AnalysisTaskFMD::CreateOutPutData() \n");
  
  fOutputList = new TList();//(TList*)GetOutputData(0);
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  
  fArray     = new TObjArray();
  fArray->SetName("FMD");
  fArray->SetOwner();
  TH1F* hEdist = 0;
 
  for(Int_t nEta = 0; nEta <= pars->GetNetaBins()+1; nEta++) {
    TObjArray* etaArray = new TObjArray();
    fArray->AddAtAndExpand(etaArray,nEta);
    for(Int_t det =1; det<=3;det++)
      {
	TObjArray* detArray = new TObjArray();
	detArray->SetName(Form("FMD%d",det));
	etaArray->AddAtAndExpand(detArray,det);
	Int_t nRings = (det==1 ? 1 : 2);
	for(Int_t ring = 0;ring<nRings;ring++)
	  {
	    Char_t ringChar = (ring == 0 ? 'I' : 'O');
	    hEdist = new TH1F(Form("FMD%d%c_etabin%d",det,ringChar,nEta),Form("FMD%d%c_etabin%d",det,ringChar,nEta),200,0,6);
	    hEdist->SetXTitle("#Delta E / E_{MIP}");
	    fOutputList->Add(hEdist);
	    detArray->AddAtAndExpand(hEdist,ring);
	  } 
      }
    
  }
  
  fZvtxDist  = new TH1F("ZvtxDist","Vertex distribution",100,-30,30);
  fZvtxDist->SetXTitle("z vertex");
  fOutputList->Add(fZvtxDist);
}

//____________________________________________________________________
void AliFMDAnalysisTaskCollector::UserExec(Option_t */*option*/)
{
  
  AliESDEvent* esd = (AliESDEvent*)InputEvent();
  AliESD* old = esd->GetAliESDOld();
  if (old) {
    esd->CopyFromOldESD();
  }
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  
  Bool_t trigger = pars->IsEventTriggered(esd);
  if(!trigger)
    return;
  Double_t vertex[3];
  
  pars->GetVertex(esd,vertex);
  if(vertex[0] == 0 && vertex[1] == 0 && vertex[2] == 0)
    return;
  
  fZvtxDist->Fill(vertex[2]);
  
  if(TMath::Abs(vertex[2]) > pars->GetVtxCutZ())
    return;
  
  AliESDFMD* fmd = esd->GetFMDData();
  if (!fmd) return;
  
  for(UShort_t det=1;det<=3;det++) {
      
    Int_t nRings = (det==1 ? 1 : 2);
    for (UShort_t ir = 0; ir < nRings; ir++) {
  
      Char_t   ring = (ir == 0 ? 'I' : 'O');
      UShort_t nsec = (ir == 0 ? 20  : 40);
      UShort_t nstr = (ir == 0 ? 512 : 256);
      TH2F* hBg = pars->GetBackgroundCorrection(det,ring,0);
      
      for(UShort_t sec =0; sec < nsec;  sec++)  {
	for(UShort_t strip = 0; strip < nstr; strip++) {
	  
	  
	  Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
	  if(mult == AliESDFMD::kInvalidMult || mult == 0) continue;
	  
	  Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]);
	  
	  Int_t nEta = hBg->GetXaxis()->FindBin(eta);
	  
	  TObjArray* etaArray = (TObjArray*)fArray->At(nEta);
	  TObjArray* detArray = (TObjArray*)etaArray->At(det);
	  TH1F* Edist = (TH1F*)detArray->At(ir);
	  
	  Edist->Fill(mult);
	  
	}
      }
    }
  }
  
  PostData(1, fOutputList); 
  
}

//____________________________________________________________________
void AliFMDAnalysisTaskCollector::ReadFromFile(const Char_t* filename, Bool_t store, Int_t speciesOption)  {

  //speciesOption:
  //0: p+p Landau fit
  //1: Pb+Pb triple landau convolution fit
  
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  pars->Init(kTRUE,AliFMDAnaParameters::kBackgroundCorrection);
  
  TFile fin(filename,"UPDATE");
  
  TList* list = (TList*)fin.Get("energyDist");
  
  AliFMDAnaCalibEnergyDistribution* EnergyDist = new AliFMDAnaCalibEnergyDistribution();
  
  EnergyDist->SetNetaBins(pars->GetNetaBins());
  EnergyDist->SetEtaLimits(pars->GetEtaMin(),pars->GetEtaMax());
    
  for(Int_t nEta = 1; nEta <= pars->GetNetaBins(); nEta++) {
  
    for(Int_t det = 1; det<=3; det++) {
      Int_t nRings  =  (det == 1 ? 1 : 2);
      for(Int_t ring = 0;ring<nRings; ring++) {
	Char_t ringChar = (ring == 0 ? 'I' : 'O');
	
	TH1F* hEdist = (TH1F*)list->FindObject(Form("FMD%d%c_etabin%d",det,ringChar,nEta));
	TF1* fitFunc = 0 ;
	
	if(hEdist->GetEntries() != 0) {
	  
	  hEdist->GetXaxis()->SetRangeUser(0.2,hEdist->GetXaxis()->GetXmax());
	  
	  if(speciesOption == 0)
	    fitFunc =  new TF1("FMDfitFunc","landau",hEdist->GetBinCenter(hEdist->GetMaximumBin())-0.2,3);
	  if(speciesOption == 1) {
	    fitFunc = new TF1("FMDfitFunc",TripleLandau,hEdist->GetBinCenter(hEdist->GetMaximumBin())-0.2,5,5);
	    fitFunc->SetParNames("constant","MPV","sigma","2-Mip weight","3-Mip weight");
	    fitFunc->SetParameters(10,0.8,0.1,0.05,0.01);
	    fitFunc->SetParLimits(1,0.6,1.2);
	    fitFunc->SetParLimits(3,0,1);
	    fitFunc->SetParLimits(4,0,1);
	    
	  }
	  
	  
	  hEdist->Fit(fitFunc,"","",hEdist->GetBinCenter(hEdist->GetMaximumBin())-0.2,3);
	  fitFunc->Write(Form("FMD%d%c_etabin%d_fitfunc",det,ringChar,nEta),TObject::kWriteDelete);
	  
	}
	
	TH2F* hBg = pars->GetBackgroundCorrection(det,ringChar,0);
	EnergyDist->SetEnergyDistribution(det,ringChar,hBg->GetXaxis()->GetBinCenter(nEta),hEdist);
      }
      
    }
    
  }
  
  fin.Write();
  fin.Close();
  
  if(store) {
    TFile fcalib(pars->GetPath(pars->GetEdistID() ),"RECREATE");
    EnergyDist->Write(AliFMDAnaParameters::GetEdistID());
    fcalib.Close();
  }


}

//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
c48a797f	1
	2	#include <TROOT.h>
	3	#include <TSystem.h>
	4	#include <TInterpreter.h>
	5	#include <TChain.h>
	6	#include <TFile.h>
	7	#include <TList.h>
9f55be54	8	#include <TMath.h>
c48a797f	9	#include <iostream>
	10
	11	#include "AliFMDAnalysisTaskCollector.h"
	12	#include "AliAnalysisManager.h"
	13	#include "AliESDFMD.h"
	14	#include "AliESDEvent.h"
	15	#include "AliAODEvent.h"
	16	#include "AliAODHandler.h"
	17	#include "AliMCEventHandler.h"
	18	#include "AliStack.h"
	19	#include "AliESDVertex.h"
	20	#include "AliFMDAnaParameters.h"
78f6f750	21	//#include "AliFMDGeometry.h"
9f55be54	22
9f55be54	23
c48a797f	24	ClassImp(AliFMDAnalysisTaskCollector)
c48a797f	25
9f55be54	26	//____________________________________________________________________
	27	Double_t AliFMDAnalysisTaskCollector::TripleLandau(Double_t x, Double_t par) {
	28
	29	Double_t energy = x[0];
	30	Double_t constant = par[0];
	31	Double_t mpv = par[1];
	32	Double_t sigma = par[2];
	33	Double_t alpha = par[3];
	34	Double_t beta = par[4];
	35
	36	Double_t f = constant*(TMath::Landau(energy,mpv,sigma,kTRUE)+
	37	alphaTMath::Landau(energy,2mpv+2sigmaTMath::Log(2),2*sigma,kTRUE)+
	38	betaTMath::Landau(energy,3mpv+3sigmaTMath::Log(3),3*sigma,kTRUE) );
	39
	40	return f;
	41	}
	42	//____________________________________________________________________
c48a797f	43
	44	AliFMDAnalysisTaskCollector::AliFMDAnalysisTaskCollector()
	45	: fDebug(0),
c48a797f	46	fOutputList(0),
c48a797f	47	fArray(0),
c48a797f	48	fZvtxDist(0)
	49	{
	50	// Default constructor
9f55be54	51
9f55be54	52
c48a797f	53	}
	54	//____________________________________________________________________
	55	AliFMDAnalysisTaskCollector::AliFMDAnalysisTaskCollector(const char* name):
9f55be54	56	AliAnalysisTaskSE(name),
c48a797f	57	fDebug(0),
c48a797f	58	fOutputList(0),
c48a797f	59	fArray(0),
c48a797f	60	fZvtxDist(0)
	61	{
	62	// Default constructor
9f55be54	63
9f55be54	64	DefineOutput(1, TList::Class());
c48a797f	65	}
c48a797f	66	//____________________________________________________________________
9f55be54	67	void AliFMDAnalysisTaskCollector::UserCreateOutputObjects()
c48a797f	68	{
	69	// Create the output container
	70	printf("AnalysisTaskFMD::CreateOutPutData() \n");
	71
	72	fOutputList = new TList();//(TList*)GetOutputData(0);
6289b3e8	73	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
c48a797f	74
	75	fArray = new TObjArray();
	76	fArray->SetName("FMD");
	77	fArray->SetOwner();
	78	TH1F* hEdist = 0;
6289b3e8	79
	80	for(Int_t nEta = 0; nEta <= pars->GetNetaBins()+1; nEta++) {
	81	TObjArray* etaArray = new TObjArray();
	82	fArray->AddAtAndExpand(etaArray,nEta);
	83	for(Int_t det =1; det<=3;det++)
	84	{
	85	TObjArray* detArray = new TObjArray();
	86	detArray->SetName(Form("FMD%d",det));
	87	etaArray->AddAtAndExpand(detArray,det);
	88	Int_t nRings = (det==1 ? 1 : 2);
	89	for(Int_t ring = 0;ring<nRings;ring++)
	90	{
	91	Char_t ringChar = (ring == 0 ? 'I' : 'O');
	92	hEdist = new TH1F(Form("FMD%d%c_etabin%d",det,ringChar,nEta),Form("FMD%d%c_etabin%d",det,ringChar,nEta),200,0,6);
	93	hEdist->SetXTitle("#Delta E / E_{MIP}");
	94	fOutputList->Add(hEdist);
	95	detArray->AddAtAndExpand(hEdist,ring);
	96	}
	97	}
	98
	99	}
c48a797f	100
	101	fZvtxDist = new TH1F("ZvtxDist","Vertex distribution",100,-30,30);
	102	fZvtxDist->SetXTitle("z vertex");
c48a797f	103	fOutputList->Add(fZvtxDist);
c48a797f	104	}
c48a797f	105
c48a797f	106	//____________________________________________________________________
9f55be54	107	void AliFMDAnalysisTaskCollector::UserExec(Option_t /option*/)
c48a797f	108	{
c48a797f	109
9f55be54	110	AliESDEvent* esd = (AliESDEvent*)InputEvent();
9f55be54	111	AliESD* old = esd->GetAliESDOld();
c48a797f	112	if (old) {
9f55be54	113	esd->CopyFromOldESD();
c48a797f	114	}
bb8a464f	115	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
c48a797f	116
9f55be54	117	Bool_t trigger = pars->IsEventTriggered(esd);
	118	if(!trigger)
	119	return;
c48a797f	120	Double_t vertex[3];
bb8a464f	121
9f55be54	122	pars->GetVertex(esd,vertex);
	123	if(vertex[0] == 0 && vertex[1] == 0 && vertex[2] == 0)
	124	return;
	125
c48a797f	126	fZvtxDist->Fill(vertex[2]);
bb8a464f	127
	128	if(TMath::Abs(vertex[2]) > pars->GetVtxCutZ())
	129	return;
	130
9f55be54	131	AliESDFMD* fmd = esd->GetFMDData();
c48a797f	132	if (!fmd) return;
	133
	134	for(UShort_t det=1;det<=3;det++) {
9f55be54	135
c48a797f	136	Int_t nRings = (det==1 ? 1 : 2);
c48a797f	137	for (UShort_t ir = 0; ir < nRings; ir++) {
6289b3e8	138
c48a797f	139	Char_t ring = (ir == 0 ? 'I' : 'O');
	140	UShort_t nsec = (ir == 0 ? 20 : 40);
	141	UShort_t nstr = (ir == 0 ? 512 : 256);
6289b3e8	142	TH2F* hBg = pars->GetBackgroundCorrection(det,ring,0);
6289b3e8	143
c48a797f	144	for(UShort_t sec =0; sec < nsec; sec++) {
c48a797f	145	for(UShort_t strip = 0; strip < nstr; strip++) {
cc066cb9	146
	147
	148	Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
	149	if(mult == AliESDFMD::kInvalidMult \|\| mult == 0) continue;
cc066cb9	150
78f6f750	151	Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]);
9f55be54	152
6289b3e8	153	Int_t nEta = hBg->GetXaxis()->FindBin(eta);
cc066cb9	154
6289b3e8	155	TObjArray* etaArray = (TObjArray*)fArray->At(nEta);
	156	TObjArray* detArray = (TObjArray*)etaArray->At(det);
	157	TH1F* Edist = (TH1F*)detArray->At(ir);
cc066cb9	158
c48a797f	159	Edist->Fill(mult);
	160
	161	}
	162	}
	163	}
	164	}
	165
9f55be54	166	PostData(1, fOutputList);
c48a797f	167
c48a797f	168	}
9f55be54	169
c48a797f	170	//____________________________________________________________________
9f55be54	171	void AliFMDAnalysisTaskCollector::ReadFromFile(const Char_t* filename, Bool_t store, Int_t speciesOption) {
	172
	173	//speciesOption:
	174	//0: p+p Landau fit
	175	//1: Pb+Pb triple landau convolution fit
	176
	177	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
	178	pars->Init(kTRUE,AliFMDAnaParameters::kBackgroundCorrection);
	179
	180	TFile fin(filename,"UPDATE");
	181
	182	TList* list = (TList*)fin.Get("energyDist");
	183
	184	AliFMDAnaCalibEnergyDistribution* EnergyDist = new AliFMDAnaCalibEnergyDistribution();
	185
	186	EnergyDist->SetNetaBins(pars->GetNetaBins());
	187	EnergyDist->SetEtaLimits(pars->GetEtaMin(),pars->GetEtaMax());
	188
	189	for(Int_t nEta = 1; nEta <= pars->GetNetaBins(); nEta++) {
	190
	191	for(Int_t det = 1; det<=3; det++) {
	192	Int_t nRings = (det == 1 ? 1 : 2);
	193	for(Int_t ring = 0;ring<nRings; ring++) {
	194	Char_t ringChar = (ring == 0 ? 'I' : 'O');
	195
	196	TH1F* hEdist = (TH1F*)list->FindObject(Form("FMD%d%c_etabin%d",det,ringChar,nEta));
	197	TF1* fitFunc = 0 ;
	198
	199	if(hEdist->GetEntries() != 0) {
	200
	201	hEdist->GetXaxis()->SetRangeUser(0.2,hEdist->GetXaxis()->GetXmax());
	202
	203	if(speciesOption == 0)
	204	fitFunc = new TF1("FMDfitFunc","landau",hEdist->GetBinCenter(hEdist->GetMaximumBin())-0.2,3);
	205	if(speciesOption == 1) {
	206	fitFunc = new TF1("FMDfitFunc",TripleLandau,hEdist->GetBinCenter(hEdist->GetMaximumBin())-0.2,5,5);
	207	fitFunc->SetParNames("constant","MPV","sigma","2-Mip weight","3-Mip weight");
	208	fitFunc->SetParameters(10,0.8,0.1,0.05,0.01);
	209	fitFunc->SetParLimits(1,0.6,1.2);
	210	fitFunc->SetParLimits(3,0,1);
	211	fitFunc->SetParLimits(4,0,1);
	212
	213	}
	214
	215
	216	hEdist->Fit(fitFunc,"","",hEdist->GetBinCenter(hEdist->GetMaximumBin())-0.2,3);
	217	fitFunc->Write(Form("FMD%d%c_etabin%d_fitfunc",det,ringChar,nEta),TObject::kWriteDelete);
	218
	219	}
	220
	221	TH2F* hBg = pars->GetBackgroundCorrection(det,ringChar,0);
	222	EnergyDist->SetEnergyDistribution(det,ringChar,hBg->GetXaxis()->GetBinCenter(nEta),hEdist);
	223	}
	224
c48a797f	225	}
bb8a464f	226
bb8a464f	227	}
bb8a464f	228
9f55be54	229	fin.Write();
	230	fin.Close();
	231
	232	if(store) {
0b0a4ae5	233	TFile fcalib(pars->GetPath(pars->GetEdistID() ),"RECREATE");
9f55be54	234	EnergyDist->Write(AliFMDAnaParameters::GetEdistID());
	235	fcalib.Close();
	236	}
	237
	238
	239
bb8a464f	240	}
9f55be54	241
c48a797f	242	//____________________________________________________________________
	243	//
	244	// EOF
	245	//