[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 <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)


AliFMDAnalysisTaskCollector::AliFMDAnalysisTaskCollector()
: fDebug(0),
  fChain(0x0),
  fESD(0x0),
  fOutputList(0),
  fArray(0),
  fEdistHist(0),
  fZvtxDist(0)
{
  // Default constructor
  DefineInput (0, TChain::Class());
  DefineOutput(0, TList::Class());
}
//____________________________________________________________________
AliFMDAnalysisTaskCollector::AliFMDAnalysisTaskCollector(const char* name):
    AliAnalysisTask(name, "AnalysisTaskFMD"),
    fDebug(0),
    fChain(0x0),
    fESD(0x0),
    fOutputList(0),
    fArray(0),
    fEdistHist(0),
    fZvtxDist(0)
{
  // Default constructor
  DefineInput (0, TChain::Class());
  DefineOutput(0, TList::Class());
}
//____________________________________________________________________
void AliFMDAnalysisTaskCollector::CreateOutputObjects()
{
  // 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(fArray);
  fOutputList->Add(fZvtxDist);
}
//____________________________________________________________________
void AliFMDAnalysisTaskCollector::Init()
{
  // Initialization
  printf("AnalysisTaskFMD::Init() \n");
  

}
//____________________________________________________________________
void AliFMDAnalysisTaskCollector::ConnectInputData(Option_t */*option*/)
{
  fChain = (TChain*)GetInputData(0);
  fESD = new AliESDEvent();
  fESD->ReadFromTree(fChain);
  
}
//____________________________________________________________________
void AliFMDAnalysisTaskCollector::Exec(Option_t */*option*/)
{
  AliESD* old = fESD->GetAliESDOld();
  if (old) {
    fESD->CopyFromOldESD();
  }
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  
  Double_t vertex[3];
  
  GetVertex(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 = fESD->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;
	  AliFMDGeometry* geo = AliFMDGeometry::Instance();
  
	  Double_t x,y,z;
	  geo->Detector2XYZ(det,ring,sec,strip,x,y,z);
	  
	  Double_t r = TMath::Sqrt(x*x+y*y);
	  
	  Double_t z_real      = z-vertex[2];
	  Double_t theta       = TMath::ATan2(r,z_real);
	  // std::cout<<"From EtaFromStrip "<<theta<<std::endl;
	  Double_t eta         =  -1*TMath::Log(TMath::Tan(0.5*theta));
	  //Float_t eta = fmd->Eta(det,ring,sec,strip);
	  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(0, fOutputList); 
  
}
//____________________________________________________________________
void AliFMDAnalysisTaskCollector::Terminate(Option_t */*option*/)
{
  /*
  for(UShort_t det=1;det<=3;det++) {
    TObjArray* detArray = (TObjArray*)fArray->At(det);
    Int_t nRings = (det==1 ? 1 : 2);
    for (UShort_t ir = 0; ir < nRings; ir++) {
      TH1F* hEdist = (TH1F*)detArray->At(ir);
      hEdist->SetAxisRange(0.4,hEdist->GetXaxis()->GetXmax());
      Float_t max = hEdist->GetBinCenter(hEdist->GetMaximumBin());
      hEdist->Fit("landau","","",max-0.1,2*max);
    }
  }
  */
}
//_____________________________________________________________________
void AliFMDAnalysisTaskCollector::GetVertex(Double_t* vertexXYZ) 
{
  const AliESDVertex* vertex = 0;
  vertex = fESD->GetPrimaryVertex();
  if (!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))        
    vertex = fESD->GetPrimaryVertexSPD();
  if (!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))        
    vertex = fESD->GetPrimaryVertexTPC();
  
  if (!vertex || (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))    
    vertex = fESD->GetVertex();
  
  if (vertex && (vertexXYZ[0] != 0 || vertexXYZ[1] != 0 || vertexXYZ[2] != 0)) {
    vertex->GetXYZ(vertexXYZ);
    //std::cout<<vertex->GetName()<<"   "<< vertex->GetTitle() <<"   "<< vertex->GetZv()<<std::endl;
    return;
  }
  else if (fESD->GetESDTZERO()) { 
    vertexXYZ[0] = 0;
    vertexXYZ[1] = 0;
    vertexXYZ[2] = fESD->GetT0zVertex();
    
    return;
  }
  
  return;
}
//____________________________________________________________________
//
// 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>
	8	#include <iostream>
	9
	10	#include "AliFMDAnalysisTaskCollector.h"
	11	#include "AliAnalysisManager.h"
	12	#include "AliESDFMD.h"
	13	#include "AliESDEvent.h"
	14	#include "AliAODEvent.h"
	15	#include "AliAODHandler.h"
	16	#include "AliMCEventHandler.h"
	17	#include "AliStack.h"
	18	#include "AliESDVertex.h"
	19	#include "AliFMDAnaParameters.h"
cc066cb9	20	#include "AliFMDGeometry.h"
c48a797f	21	ClassImp(AliFMDAnalysisTaskCollector)
	22
	23
	24	AliFMDAnalysisTaskCollector::AliFMDAnalysisTaskCollector()
	25	: fDebug(0),
	26	fChain(0x0),
	27	fESD(0x0),
c48a797f	28	fOutputList(0),
	29	fArray(0),
	30	fEdistHist(0),
	31	fZvtxDist(0)
	32	{
	33	// Default constructor
	34	DefineInput (0, TChain::Class());
	35	DefineOutput(0, TList::Class());
	36	}
	37	//____________________________________________________________________
	38	AliFMDAnalysisTaskCollector::AliFMDAnalysisTaskCollector(const char* name):
	39	AliAnalysisTask(name, "AnalysisTaskFMD"),
	40	fDebug(0),
	41	fChain(0x0),
	42	fESD(0x0),
c48a797f	43	fOutputList(0),
	44	fArray(0),
	45	fEdistHist(0),
	46	fZvtxDist(0)
	47	{
	48	// Default constructor
	49	DefineInput (0, TChain::Class());
	50	DefineOutput(0, TList::Class());
	51	}
	52	//____________________________________________________________________
	53	void AliFMDAnalysisTaskCollector::CreateOutputObjects()
	54	{
	55	// Create the output container
	56	printf("AnalysisTaskFMD::CreateOutPutData() \n");
	57
	58	fOutputList = new TList();//(TList*)GetOutputData(0);
6289b3e8	59	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
c48a797f	60
	61	fArray = new TObjArray();
	62	fArray->SetName("FMD");
	63	fArray->SetOwner();
	64	TH1F* hEdist = 0;
6289b3e8	65
	66	for(Int_t nEta = 0; nEta <= pars->GetNetaBins()+1; nEta++) {
	67	TObjArray* etaArray = new TObjArray();
	68	fArray->AddAtAndExpand(etaArray,nEta);
	69	for(Int_t det =1; det<=3;det++)
	70	{
	71	TObjArray* detArray = new TObjArray();
	72	detArray->SetName(Form("FMD%d",det));
	73	etaArray->AddAtAndExpand(detArray,det);
	74	Int_t nRings = (det==1 ? 1 : 2);
	75	for(Int_t ring = 0;ring<nRings;ring++)
	76	{
	77	Char_t ringChar = (ring == 0 ? 'I' : 'O');
	78	hEdist = new TH1F(Form("FMD%d%c_etabin%d",det,ringChar,nEta),Form("FMD%d%c_etabin%d",det,ringChar,nEta),200,0,6);
	79	hEdist->SetXTitle("#Delta E / E_{MIP}");
	80	fOutputList->Add(hEdist);
	81	detArray->AddAtAndExpand(hEdist,ring);
	82	}
	83	}
	84
	85	}
c48a797f	86
	87	fZvtxDist = new TH1F("ZvtxDist","Vertex distribution",100,-30,30);
	88	fZvtxDist->SetXTitle("z vertex");
b82e76e0	89	//fOutputList->Add(fArray);
c48a797f	90	fOutputList->Add(fZvtxDist);
	91	}
	92	//____________________________________________________________________
	93	void AliFMDAnalysisTaskCollector::Init()
	94	{
	95	// Initialization
	96	printf("AnalysisTaskFMD::Init() \n");
	97
	98
	99	}
	100	//____________________________________________________________________
	101	void AliFMDAnalysisTaskCollector::ConnectInputData(Option_t /option*/)
	102	{
	103	fChain = (TChain*)GetInputData(0);
	104	fESD = new AliESDEvent();
	105	fESD->ReadFromTree(fChain);
	106
	107	}
	108	//____________________________________________________________________
	109	void AliFMDAnalysisTaskCollector::Exec(Option_t /option*/)
	110	{
	111	AliESD* old = fESD->GetAliESDOld();
	112	if (old) {
	113	fESD->CopyFromOldESD();
	114	}
bb8a464f	115	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
c48a797f	116
c48a797f	117	Double_t vertex[3];
bb8a464f	118
bb8a464f	119	GetVertex(vertex);
6289b3e8	120	//if(vertex[0] == 0 && vertex[1] == 0 && vertex[2] == 0)
6289b3e8	121	// return;
c48a797f	122	fZvtxDist->Fill(vertex[2]);
bb8a464f	123
	124	if(TMath::Abs(vertex[2]) > pars->GetVtxCutZ())
	125	return;
	126
c48a797f	127	AliESDFMD* fmd = fESD->GetFMDData();
	128	if (!fmd) return;
	129
6289b3e8	130
6289b3e8	131
c48a797f	132	for(UShort_t det=1;det<=3;det++) {
6289b3e8	133
6289b3e8	134
c48a797f	135	Int_t nRings = (det==1 ? 1 : 2);
c48a797f	136	for (UShort_t ir = 0; ir < nRings; ir++) {
6289b3e8	137
c48a797f	138	Char_t ring = (ir == 0 ? 'I' : 'O');
	139	UShort_t nsec = (ir == 0 ? 20 : 40);
	140	UShort_t nstr = (ir == 0 ? 512 : 256);
6289b3e8	141	TH2F* hBg = pars->GetBackgroundCorrection(det,ring,0);
6289b3e8	142
c48a797f	143	for(UShort_t sec =0; sec < nsec; sec++) {
c48a797f	144	for(UShort_t strip = 0; strip < nstr; strip++) {
cc066cb9	145
	146
	147	Float_t mult = fmd->Multiplicity(det,ring,sec,strip);
	148	if(mult == AliESDFMD::kInvalidMult \|\| mult == 0) continue;
	149	AliFMDGeometry* geo = AliFMDGeometry::Instance();
	150
	151	Double_t x,y,z;
	152	geo->Detector2XYZ(det,ring,sec,strip,x,y,z);
	153
	154	Double_t r = TMath::Sqrt(xx+yy);
	155
	156	Double_t z_real = z-vertex[2];
	157	Double_t theta = TMath::ATan2(r,z_real);
	158	// std::cout<<"From EtaFromStrip "<<theta<<std::endl;
	159	Double_t eta = -1TMath::Log(TMath::Tan(0.5theta));
	160	//Float_t eta = fmd->Eta(det,ring,sec,strip);
6289b3e8	161	Int_t nEta = hBg->GetXaxis()->FindBin(eta);
cc066cb9	162
6289b3e8	163	TObjArray* etaArray = (TObjArray*)fArray->At(nEta);
	164	TObjArray* detArray = (TObjArray*)etaArray->At(det);
	165	TH1F* Edist = (TH1F*)detArray->At(ir);
cc066cb9	166
c48a797f	167	Edist->Fill(mult);
	168
	169	}
	170	}
	171	}
	172	}
	173
	174	PostData(0, fOutputList);
	175
	176	}
	177	//____________________________________________________________________
	178	void AliFMDAnalysisTaskCollector::Terminate(Option_t /option*/)
	179	{
b82e76e0	180	/*
c48a797f	181	for(UShort_t det=1;det<=3;det++) {
	182	TObjArray* detArray = (TObjArray*)fArray->At(det);
	183	Int_t nRings = (det==1 ? 1 : 2);
	184	for (UShort_t ir = 0; ir < nRings; ir++) {
	185	TH1F* hEdist = (TH1F*)detArray->At(ir);
	186	hEdist->SetAxisRange(0.4,hEdist->GetXaxis()->GetXmax());
	187	Float_t max = hEdist->GetBinCenter(hEdist->GetMaximumBin());
	188	hEdist->Fit("landau","","",max-0.1,2*max);
	189	}
	190	}
b82e76e0	191	*/
c48a797f	192	}
bb8a464f	193	//_____________________________________________________________________
	194	void AliFMDAnalysisTaskCollector::GetVertex(Double_t* vertexXYZ)
	195	{
	196	const AliESDVertex* vertex = 0;
	197	vertex = fESD->GetPrimaryVertex();
	198	if (!vertex \|\| (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))
	199	vertex = fESD->GetPrimaryVertexSPD();
	200	if (!vertex \|\| (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))
	201	vertex = fESD->GetPrimaryVertexTPC();
	202
	203	if (!vertex \|\| (vertexXYZ[0] == 0 && vertexXYZ[1] == 0 && vertexXYZ[2] == 0))
	204	vertex = fESD->GetVertex();
	205
	206	if (vertex && (vertexXYZ[0] != 0 \|\| vertexXYZ[1] != 0 \|\| vertexXYZ[2] != 0)) {
	207	vertex->GetXYZ(vertexXYZ);
	208	//std::cout<<vertex->GetName()<<" "<< vertex->GetTitle() <<" "<< vertex->GetZv()<<std::endl;
	209	return;
	210	}
	211	else if (fESD->GetESDTZERO()) {
	212	vertexXYZ[0] = 0;
	213	vertexXYZ[1] = 0;
	214	vertexXYZ[2] = fESD->GetT0zVertex();
	215
	216	return;
	217	}
	218
	219	return;
	220	}
c48a797f	221	//____________________________________________________________________
	222	//
	223	// EOF
	224	//