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

 
#include <TROOT.h>
#include <TSystem.h>
#include <TInterpreter.h>
#include <TChain.h>
#include <TFile.h>
#include <TList.h>
#include <iostream>
#include "TH1F.h"
#include "TH2F.h"
#include "AliFMDAnalysisTaskDndeta.h"
#include "AliAnalysisManager.h"
#include "AliESDFMD.h"
#include "AliESDEvent.h"
#include "AliAODEvent.h"
#include "AliAODHandler.h"
#include "AliMCEventHandler.h"
#include "AliStack.h"
#include "AliLog.h"
#include "AliESDVertex.h"
#include "TMath.h"
#include "AliFMDAnaParameters.h"
#include "AliFMDGeometry.h"
#include "AliGenEventHeader.h"
#include "TDatabasePDG.h"
#include "TParticlePDG.h"
ClassImp(AliFMDAnalysisTaskDndeta)


AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta()
: fDebug(0),
  fOutputList(0),
  fInputList(0),
  fArray(0),
  fInputArray(0),
  fVertexString(0x0),
  fNevents(),
  fNMCevents(),
  fStandalone(kTRUE),
  fMCevent(0)
{
  // Default constructor
  DefineInput (0, TList::Class());
  DefineOutput(0, TList::Class());
}
//_____________________________________________________________________
AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta(const char* name, Bool_t SE):
    AliAnalysisTask(name, "Density"),
    fDebug(0),
    fOutputList(0),
    fInputList(0),
    fArray(),
    fInputArray(0),
    fVertexString(0x0),
    fNevents(),
    fNMCevents(),
    fStandalone(kTRUE),
    fMCevent(0)
{
  fStandalone = SE;
  if(fStandalone) {
    DefineInput (0, TList::Class());
    DefineInput(1, TObjString::Class());
    DefineOutput(0, TList::Class());
    
  }
}
//_____________________________________________________________________
void AliFMDAnalysisTaskDndeta::CreateOutputObjects()
{
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  
  fArray.SetName("FMD");
  fArray.SetOwner();
  
  if(!fOutputList)
    fOutputList = new TList();
  fOutputList->SetName("BackgroundCorrected");
  
  
  TH2F* hMult = 0;
  TH1F* hPrimVertexBin = 0;
  
  
  TH2F* hBg = pars->GetBackgroundCorrection(1, 'I', 0);
  TH1F* hPrimary = new TH1F("hMultvsEta","hMultvsEta",
			    hBg->GetNbinsX(),
			    hBg->GetXaxis()->GetXmin(),
			    hBg->GetXaxis()->GetXmax());
  hPrimary->Sumw2();
  fOutputList->Add(hPrimary);
  Int_t nVtxbins = pars->GetNvtxBins();
  
  
  for(Int_t det =1; det<=3;det++)
    {
      TObjArray* detArray = new TObjArray();
      detArray->SetName(Form("FMD%d",det));
      fArray.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');
	  Int_t  nSec     = (ring == 0 ? 20 : 40);
	  
	  TObjArray* vtxArray = new TObjArray();
	  vtxArray->SetName(Form("FMD%d%c",det,ringChar));
	  detArray->AddAtAndExpand(vtxArray,ring);
	  for(Int_t i = 0; i< nVtxbins; i++) {
	    TH2F* hBg = pars->GetBackgroundCorrection(det, ringChar, i);
	    hMult  = new TH2F(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),
			      hBg->GetNbinsX(),
			      hBg->GetXaxis()->GetXmin(),
			      hBg->GetXaxis()->GetXmax(),
			      nSec, 0, 2*TMath::Pi());
	    hMult->Sumw2();
	    fOutputList->Add(hMult);
	    vtxArray->AddAtAndExpand(hMult,i);
	    
	  }
	} 
    }
  
  for(Int_t i = 0; i< nVtxbins; i++) {
   
    hPrimVertexBin = new TH1F(Form("primmult_vtxbin%d",i),
			      Form("primmult_vtxbin%d",i),
			      hBg->GetNbinsX(),
			      hBg->GetXaxis()->GetXmin(),
			      hBg->GetXaxis()->GetXmax());
    hPrimVertexBin->Sumw2();
    fOutputList->Add(hPrimVertexBin);
    
  }
  
  fNevents.SetBins(nVtxbins,0,nVtxbins);
  fNevents.SetName("nEvents");
  fNMCevents.SetBins(nVtxbins,0,nVtxbins);
  fNMCevents.SetName("nMCEvents");
  fOutputList->Add(&fNevents);
  fOutputList->Add(&fNMCevents);
  
}
//_____________________________________________________________________
void AliFMDAnalysisTaskDndeta::ConnectInputData(Option_t */*option*/)
{
  if(fStandalone) {
    fInputList   = (TList*)GetInputData(0);
    fVertexString = (TObjString*)GetInputData(1);
  }
}
//_____________________________________________________________________
void AliFMDAnalysisTaskDndeta::Exec(Option_t */*option*/)
{
  Int_t vtxbin   = fVertexString->GetString().Atoi();
  fNevents.Fill(vtxbin);
  for(UShort_t det=1;det<=3;det++) {
    //TObjArray* detInputArray = (TObjArray*)fInputArray->At(det);
    TObjArray* detArray = (TObjArray*)fArray.At(det);
    Int_t nRings = (det==1 ? 1 : 2);
    for (UShort_t ir = 0; ir < nRings; ir++) {
      Char_t ringChar = (ir == 0 ? 'I' : 'O');
      //TObjArray* vtxInputArray = (TObjArray*)detInputArray->At(ir);
      TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
      TH2F* hMultTotal = (TH2F*)vtxArray->At(vtxbin);
      
      
      TH2F* hMultInput = (TH2F*)fInputList->FindObject(Form("mult_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
      
      hMultTotal->Add(hMultInput);
      
      
    }
  }
  
  if(fMCevent)
    ProcessPrimary();
  
  if(fStandalone) {
    PostData(0, fOutputList); 
  }
  
}
//_____________________________________________________________________
void AliFMDAnalysisTaskDndeta::Terminate(Option_t */*option*/) {
  
  
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  
  Int_t nVtxbins = pars->GetNvtxBins();
  
  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++) {
      TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
      Char_t ringChar = (ir == 0 ? 'I' : 'O');
      for(Int_t i =0; i<nVtxbins; i++) {
	TH2F* hMultTotal = (TH2F*)vtxArray->At(i);
	TH1D* hMultProj   = hMultTotal->ProjectionX(Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ringChar,i),1,hMultTotal->GetNbinsY());
	fOutputList->Add(hMultProj);
      }
    }
  }
  
}
//_____________________________________________________________________
void AliFMDAnalysisTaskDndeta::ProcessPrimary() {

  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  
  AliMCParticle* particle = 0;
  AliStack* stack = fMCevent->Stack();
  
  TH1F* hPrimary = (TH1F*)fOutputList->FindObject("hMultvsEta");
  
  Bool_t firstTrack = kTRUE;
  Int_t nTracks = fMCevent->GetNumberOfTracks();
  for(Int_t i = 0 ;i<nTracks;i++) {
    particle = fMCevent->GetTrack(i);
    if(!particle)
      continue;
    if(TMath::Abs(particle->Zv()) > pars->GetVtxCutZ())
      continue;
    
    if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
      hPrimary->Fill(particle->Eta());
      Double_t delta           = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
      Double_t vertexBinDouble = (particle->Zv() + pars->GetVtxCutZ()) / delta;
      Int_t    vertexBin       = (Int_t)vertexBinDouble;

      TH1F* hPrimVtxBin = (TH1F*)fOutputList->FindObject(Form("primmult_vtxbin%d",vertexBin));
      hPrimVtxBin->Fill(particle->Eta());
      if(firstTrack) {
	fNMCevents.Fill(vertexBin);
	firstTrack = kFALSE;
      }
    }
      
  }
  
  
}
//_____________________________________________________________________
//
//
// EOF
Commit	Line	Data
7c3e5162	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>
bb8a464f	9	#include "TH1F.h"
7c3e5162	10	#include "TH2F.h"
	11	#include "AliFMDAnalysisTaskDndeta.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 "AliLog.h"
	20	#include "AliESDVertex.h"
	21	#include "TMath.h"
	22	#include "AliFMDAnaParameters.h"
	23	#include "AliFMDGeometry.h"
bb8a464f	24	#include "AliGenEventHeader.h"
6289b3e8	25	#include "TDatabasePDG.h"
6289b3e8	26	#include "TParticlePDG.h"
7c3e5162	27	ClassImp(AliFMDAnalysisTaskDndeta)
	28
	29
	30	AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta()
	31	: fDebug(0),
	32	fOutputList(0),
	33	fInputList(0),
	34	fArray(0),
	35	fInputArray(0),
	36	fVertexString(0x0),
	37	fNevents(),
bb8a464f	38	fNMCevents(),
	39	fStandalone(kTRUE),
	40	fMCevent(0)
7c3e5162	41	{
	42	// Default constructor
	43	DefineInput (0, TList::Class());
	44	DefineOutput(0, TList::Class());
	45	}
	46	//_____________________________________________________________________
	47	AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta(const char* name, Bool_t SE):
	48	AliAnalysisTask(name, "Density"),
	49	fDebug(0),
	50	fOutputList(0),
	51	fInputList(0),
	52	fArray(),
	53	fInputArray(0),
	54	fVertexString(0x0),
	55	fNevents(),
bb8a464f	56	fNMCevents(),
	57	fStandalone(kTRUE),
	58	fMCevent(0)
7c3e5162	59	{
	60	fStandalone = SE;
	61	if(fStandalone) {
	62	DefineInput (0, TList::Class());
	63	DefineInput(1, TObjString::Class());
	64	DefineOutput(0, TList::Class());
	65
	66	}
	67	}
	68	//_____________________________________________________________________
	69	void AliFMDAnalysisTaskDndeta::CreateOutputObjects()
	70	{
	71	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
	72
	73	fArray.SetName("FMD");
	74	fArray.SetOwner();
	75
	76	if(!fOutputList)
	77	fOutputList = new TList();
	78	fOutputList->SetName("BackgroundCorrected");
	79
	80
	81	TH2F* hMult = 0;
bb8a464f	82	TH1F* hPrimVertexBin = 0;
7c3e5162	83
bb8a464f	84
	85	TH2F* hBg = pars->GetBackgroundCorrection(1, 'I', 0);
	86	TH1F* hPrimary = new TH1F("hMultvsEta","hMultvsEta",
	87	hBg->GetNbinsX(),
	88	hBg->GetXaxis()->GetXmin(),
	89	hBg->GetXaxis()->GetXmax());
	90	hPrimary->Sumw2();
	91	fOutputList->Add(hPrimary);
7c3e5162	92	Int_t nVtxbins = pars->GetNvtxBins();
7c3e5162	93
bb8a464f	94
7c3e5162	95	for(Int_t det =1; det<=3;det++)
	96	{
	97	TObjArray* detArray = new TObjArray();
	98	detArray->SetName(Form("FMD%d",det));
	99	fArray.AddAtAndExpand(detArray,det);
	100	Int_t nRings = (det==1 ? 1 : 2);
	101	for(Int_t ring = 0;ring<nRings;ring++)
	102	{
	103	Char_t ringChar = (ring == 0 ? 'I' : 'O');
	104	Int_t nSec = (ring == 0 ? 20 : 40);
	105
	106	TObjArray* vtxArray = new TObjArray();
	107	vtxArray->SetName(Form("FMD%d%c",det,ringChar));
	108	detArray->AddAtAndExpand(vtxArray,ring);
	109	for(Int_t i = 0; i< nVtxbins; i++) {
	110	TH2F* hBg = pars->GetBackgroundCorrection(det, ringChar, i);
	111	hMult = new TH2F(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),
	112	hBg->GetNbinsX(),
	113	hBg->GetXaxis()->GetXmin(),
	114	hBg->GetXaxis()->GetXmax(),
	115	nSec, 0, 2*TMath::Pi());
	116	hMult->Sumw2();
	117	fOutputList->Add(hMult);
	118	vtxArray->AddAtAndExpand(hMult,i);
	119
	120	}
	121	}
	122	}
	123
bb8a464f	124	for(Int_t i = 0; i< nVtxbins; i++) {
	125
	126	hPrimVertexBin = new TH1F(Form("primmult_vtxbin%d",i),
	127	Form("primmult_vtxbin%d",i),
	128	hBg->GetNbinsX(),
	129	hBg->GetXaxis()->GetXmin(),
	130	hBg->GetXaxis()->GetXmax());
	131	hPrimVertexBin->Sumw2();
	132	fOutputList->Add(hPrimVertexBin);
	133
	134	}
	135
7c3e5162	136	fNevents.SetBins(nVtxbins,0,nVtxbins);
7c3e5162	137	fNevents.SetName("nEvents");
bb8a464f	138	fNMCevents.SetBins(nVtxbins,0,nVtxbins);
bb8a464f	139	fNMCevents.SetName("nMCEvents");
7c3e5162	140	fOutputList->Add(&fNevents);
bb8a464f	141	fOutputList->Add(&fNMCevents);
7c3e5162	142
	143	}
	144	//_____________________________________________________________________
	145	void AliFMDAnalysisTaskDndeta::ConnectInputData(Option_t /option*/)
	146	{
	147	if(fStandalone) {
	148	fInputList = (TList*)GetInputData(0);
	149	fVertexString = (TObjString*)GetInputData(1);
	150	}
	151	}
	152	//_____________________________________________________________________
	153	void AliFMDAnalysisTaskDndeta::Exec(Option_t /option*/)
	154	{
	155	Int_t vtxbin = fVertexString->GetString().Atoi();
	156	fNevents.Fill(vtxbin);
7c3e5162	157	for(UShort_t det=1;det<=3;det++) {
	158	//TObjArray* detInputArray = (TObjArray*)fInputArray->At(det);
	159	TObjArray* detArray = (TObjArray*)fArray.At(det);
	160	Int_t nRings = (det==1 ? 1 : 2);
	161	for (UShort_t ir = 0; ir < nRings; ir++) {
	162	Char_t ringChar = (ir == 0 ? 'I' : 'O');
	163	//TObjArray* vtxInputArray = (TObjArray*)detInputArray->At(ir);
	164	TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
	165	TH2F* hMultTotal = (TH2F*)vtxArray->At(vtxbin);
	166
	167
	168	TH2F* hMultInput = (TH2F*)fInputList->FindObject(Form("mult_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
1282ce49	169
7c3e5162	170	hMultTotal->Add(hMultInput);
	171
	172
	173	}
	174	}
bb8a464f	175
	176	if(fMCevent)
	177	ProcessPrimary();
	178
7c3e5162	179	if(fStandalone) {
	180	PostData(0, fOutputList);
	181	}
	182
	183	}
	184	//_____________________________________________________________________
	185	void AliFMDAnalysisTaskDndeta::Terminate(Option_t /option*/) {
	186
1282ce49	187
7c3e5162	188	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
	189
	190	Int_t nVtxbins = pars->GetNvtxBins();
	191
	192	for(UShort_t det=1;det<=3;det++) {
	193	TObjArray* detArray = (TObjArray*)fArray.At(det);
	194	Int_t nRings = (det==1 ? 1 : 2);
	195	for (UShort_t ir = 0; ir < nRings; ir++) {
	196	TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
1282ce49	197	Char_t ringChar = (ir == 0 ? 'I' : 'O');
7c3e5162	198	for(Int_t i =0; i<nVtxbins; i++) {
7c3e5162	199	TH2F* hMultTotal = (TH2F*)vtxArray->At(i);
1282ce49	200	TH1D* hMultProj = hMultTotal->ProjectionX(Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ringChar,i),1,hMultTotal->GetNbinsY());
1282ce49	201	fOutputList->Add(hMultProj);
7c3e5162	202	}
	203	}
	204	}
1282ce49	205
bb8a464f	206	}
	207	//_____________________________________________________________________
	208	void AliFMDAnalysisTaskDndeta::ProcessPrimary() {
	209
	210	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
	211
	212	AliMCParticle* particle = 0;
	213	AliStack* stack = fMCevent->Stack();
	214
	215	TH1F* hPrimary = (TH1F*)fOutputList->FindObject("hMultvsEta");
	216
	217	Bool_t firstTrack = kTRUE;
	218	Int_t nTracks = fMCevent->GetNumberOfTracks();
	219	for(Int_t i = 0 ;i<nTracks;i++) {
	220	particle = fMCevent->GetTrack(i);
	221	if(!particle)
	222	continue;
	223	if(TMath::Abs(particle->Zv()) > pars->GetVtxCutZ())
	224	continue;
6289b3e8	225
bb8a464f	226	if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
	227	hPrimary->Fill(particle->Eta());
	228	Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
	229	Double_t vertexBinDouble = (particle->Zv() + pars->GetVtxCutZ()) / delta;
	230	Int_t vertexBin = (Int_t)vertexBinDouble;
6289b3e8	231
bb8a464f	232	TH1F* hPrimVtxBin = (TH1F*)fOutputList->FindObject(Form("primmult_vtxbin%d",vertexBin));
	233	hPrimVtxBin->Fill(particle->Eta());
	234	if(firstTrack) {
	235	fNMCevents.Fill(vertexBin);
	236	firstTrack = kFALSE;
	237	}
	238	}
	239
	240	}
	241
	242
	243
	244
	245
7c3e5162	246	}
	247	//_____________________________________________________________________
	248	//
	249	//
	250	// EOF