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

 
#include <TROOT.h>
#include <TSystem.h>
#include <TInterpreter.h>
#include <TChain.h>
#include <TFile.h>
#include <TList.h>
#include <iostream>
#include "TAxis.h"
#include "TH2F.h"
#include "AliFMDAnalysisTaskDensity.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 "TMath.h"
#include "AliFMDAnaParameters.h"
#include "AliFMDGeometry.h"
#include "AliFMDRing.h"

ClassImp(AliFMDAnalysisTaskDensity)

//_____________________________________________________________________
AliFMDAnalysisTaskDensity::AliFMDAnalysisTaskDensity()
: fDebug(0),
  fOutputList(),
  fArray(),
  fESD(0x0),
  fVertexString(),
  fVertex(0),
  fStandalone(kTRUE),
  fStatus(kTRUE)
{
  // Default constructor
  DefineInput (0, AliESDFMD::Class());
  DefineInput (1, AliESDVertex::Class());
  DefineOutput(0,TList::Class());
}
//_____________________________________________________________________
AliFMDAnalysisTaskDensity::AliFMDAnalysisTaskDensity(const char* name, Bool_t SE):
    AliAnalysisTask(name, "Density"),
    fDebug(0),
    fOutputList(0),
    fArray(),
    fESD(0x0),
    fVertexString(),
    fVertex(0),
    fStandalone(kTRUE),
    fStatus(kTRUE)
{
  fStandalone = SE;
  if(fStandalone) {
    DefineInput (0, AliESDFMD::Class());
    DefineInput (1, AliESDVertex::Class());
    DefineOutput(0, TList::Class());
  }
}
//_____________________________________________________________________
void AliFMDAnalysisTaskDensity::CreateOutputObjects()
{
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  
  fArray.SetName("FMD");
  fArray.SetOwner();
  if(!fOutputList)
    fOutputList = new TList();
  fOutputList->SetName("density_list");
  
  fOutputList->Add(&fArray);
  fOutputList->Add(&fVertexString);
  
  TH2F* hMult = 0;
  
  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("FMD%d%c_vtxbin%d",det,ringChar,i),Form("FMD%d%c_vtxbin%d",det,ringChar,i),
			      hBg->GetNbinsX(),
			      hBg->GetXaxis()->GetXmin(),
			      hBg->GetXaxis()->GetXmax(),
			      nSec, 0, 2*TMath::Pi());
	    
	    vtxArray->AddAtAndExpand(hMult,i);
	  }
	} 
    }
  
  
  
  
}
//_____________________________________________________________________
void AliFMDAnalysisTaskDensity::ConnectInputData(Option_t */*option*/)
{
  if(fStandalone) {
    fESD    = (AliESDFMD*)GetInputData(0);
    fVertex = (AliESDVertex*)GetInputData(1);
  }
}
//_____________________________________________________________________
void AliFMDAnalysisTaskDensity::Exec(Option_t */*option*/)
{
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  AliFMDGeometry* geo       = AliFMDGeometry::Instance();
  
  //AliESDFMD*   fmd = fESD->GetFMDData();
  
  Double_t vertex[3];
  fVertex->GetXYZ(vertex);
  // Z Vtx cut
  if( TMath::Abs(vertex[2]) > pars->GetVtxCutZ()) {
    fStatus = kFALSE;
    return;
  }
  else
    fStatus = kTRUE;
  
  Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
  Double_t vertexBinDouble = (vertex[2] + pars->GetVtxCutZ()) / delta;
  
  Int_t vtxbin = (Int_t)vertexBinDouble;

  fVertexString.SetString(Form("%d",vtxbin));
  //Reset everything
  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);
      
      TH2F* hMult   = (TH2F*)vtxArray->At(vtxbin); 
      hMult->Reset();
    }
    
  }
  
  
  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);
      
      TH2F* hMult   = (TH2F*)vtxArray->At(vtxbin);
      
      Char_t   ring = (ir == 0 ? 'I' : 'O');
      UShort_t nsec = (ir == 0 ? 20  : 40);
      UShort_t nstr = (ir == 0 ? 512 : 256);
      
      for(UShort_t sec =0; sec < nsec;  sec++)  {
	for(UShort_t strip = 0; strip < nstr; strip++) {
	  Float_t mult = fESD->Multiplicity(det,ring,sec,strip);
	  Float_t eta = fESD->Eta(det,ring,sec,strip);
	  Float_t mult_cut = 0.2;
	  if(mult == 0 || mult == AliESDFMD::kInvalidMult) continue;
	  //Particle number cut goes here...
	  Float_t nParticles = 0;
	  if(fESD->GetUniqueID() == kTRUE) {
	    //proton + proton
	    if(mult > mult_cut) 
	      nParticles = 1; 
	  }
	  else {
	    
	    //Pb+Pb
	    Float_t mpv   = pars->GetMPV(det,ring,eta);
	    Float_t sigma = pars->GetSigma(det,ring,eta);
	    Float_t alpha = pars->Get2MIPWeight(det,ring,eta);
	    Float_t beta  = pars->Get3MIPWeight(det,ring,eta);
	    
	    Float_t sumCor = TMath::Landau(mult,mpv,sigma,kTRUE)+
	      alpha*TMath::Landau(mult,2*mpv+2*sigma*TMath::Log(2),2*sigma,kTRUE)+
	      beta*TMath::Landau(mult,3*mpv+3*sigma*TMath::Log(3),3*sigma,kTRUE);
	    Float_t weight = TMath::Landau(mult,mpv,sigma,kTRUE)+
	      2*alpha*TMath::Landau(mult,2*mpv+2*sigma*TMath::Log(2),2*sigma,kTRUE)+
	      3*beta*TMath::Landau(mult,3*mpv+3*sigma*TMath::Log(3),3*sigma,kTRUE);
	    
	    
	    if(mult > 0){//mult_cut) {
	      if(sumCor) nParticles = weight / sumCor;
	      else nParticles = 1;
	    }
	    //std::cout<<sumCor<<"    "<<weight<<"    "<<"    "<<mult<<"  "<<nParticles<<std::endl;
	    
	  }
	  
	  
	  Double_t x,y,z;
	  geo->Detector2XYZ(det,ring,sec,strip,x,y,z);
	  Float_t phi = TMath::ATan2(y,x);
	  if(phi<0)
	    phi = phi+2*TMath::Pi();
	  Float_t correction = GetAcceptanceCorrection(ring,strip);
	  if(correction) nParticles = nParticles / correction;
	  hMult->Fill(eta,phi,nParticles);
	  
	  
	}
      }
      
    }
    
	
  
  }
  if(fStandalone) {
    PostData(0, fOutputList); 
  }
  
}
//_____________________________________________________________________
Float_t AliFMDAnalysisTaskDensity::GetAcceptanceCorrection(Char_t ring, UShort_t strip)
{
  AliFMDRing fmdring(ring);
  fmdring.Init();
  Float_t   rad       = fmdring.GetMaxR()-fmdring.GetMinR();
  Float_t   segment   = rad / fmdring.GetNStrips();
  Float_t   radius    = fmdring.GetMinR() + segment*strip;
  
  Float_t   basearea1 = 0.5*fmdring.GetBaseStripLength(strip)*TMath::Power(radius,2);
  Float_t   basearea2 = 0.5*fmdring.GetBaseStripLength(strip)*TMath::Power((radius-segment),2);
  Float_t   basearea  = basearea1 - basearea2;
  
  Float_t   area1     = 0.5*fmdring.GetStripLength(strip)*TMath::Power(radius,2);
  Float_t   area2     = 0.5*fmdring.GetStripLength(strip)*TMath::Power((radius-segment),2);
  Float_t   area      = area1 - area2;
  
  Float_t correction = area/basearea;
  
  return correction;
}

//
//EOF
//
Commit	Line	Data
3bb122c7	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	#include "TAxis.h"
	10	#include "TH2F.h"
	11	#include "AliFMDAnalysisTaskDensity.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 "TMath.h"
	21	#include "AliFMDAnaParameters.h"
	22	#include "AliFMDGeometry.h"
1282ce49	23	#include "AliFMDRing.h"
3bb122c7	24
	25	ClassImp(AliFMDAnalysisTaskDensity)
	26
	27	//_____________________________________________________________________
	28	AliFMDAnalysisTaskDensity::AliFMDAnalysisTaskDensity()
	29	: fDebug(0),
8dc7c4c2	30	fOutputList(),
8dc7c4c2	31	fArray(),
c78bc12b	32	fESD(0x0),
7c3e5162	33	fVertexString(),
7c3e5162	34	fVertex(0),
bb8a464f	35	fStandalone(kTRUE),
bb8a464f	36	fStatus(kTRUE)
3bb122c7	37	{
3bb122c7	38	// Default constructor
7c3e5162	39	DefineInput (0, AliESDFMD::Class());
7c3e5162	40	DefineInput (1, AliESDVertex::Class());
3bb122c7	41	DefineOutput(0,TList::Class());
	42	}
	43	//_____________________________________________________________________
7c3e5162	44	AliFMDAnalysisTaskDensity::AliFMDAnalysisTaskDensity(const char* name, Bool_t SE):
3bb122c7	45	AliAnalysisTask(name, "Density"),
3bb122c7	46	fDebug(0),
7c3e5162	47	fOutputList(0),
8dc7c4c2	48	fArray(),
c78bc12b	49	fESD(0x0),
7c3e5162	50	fVertexString(),
7c3e5162	51	fVertex(0),
bb8a464f	52	fStandalone(kTRUE),
bb8a464f	53	fStatus(kTRUE)
3bb122c7	54	{
7c3e5162	55	fStandalone = SE;
	56	if(fStandalone) {
	57	DefineInput (0, AliESDFMD::Class());
	58	DefineInput (1, AliESDVertex::Class());
	59	DefineOutput(0, TList::Class());
	60	}
3bb122c7	61	}
	62	//_____________________________________________________________________
	63	void AliFMDAnalysisTaskDensity::CreateOutputObjects()
	64	{
	65	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
3bb122c7	66
8dc7c4c2	67	fArray.SetName("FMD");
8dc7c4c2	68	fArray.SetOwner();
7c3e5162	69	if(!fOutputList)
	70	fOutputList = new TList();
	71	fOutputList->SetName("density_list");
3bb122c7	72
bb8a464f	73	fOutputList->Add(&fArray);
	74	fOutputList->Add(&fVertexString);
	75
3bb122c7	76	TH2F* hMult = 0;
	77
	78	Int_t nVtxbins = pars->GetNvtxBins();
	79
	80	for(Int_t det =1; det<=3;det++)
	81	{
	82	TObjArray* detArray = new TObjArray();
	83	detArray->SetName(Form("FMD%d",det));
8dc7c4c2	84	fArray.AddAtAndExpand(detArray,det);
3bb122c7	85	Int_t nRings = (det==1 ? 1 : 2);
	86	for(Int_t ring = 0;ring<nRings;ring++)
	87	{
	88	Char_t ringChar = (ring == 0 ? 'I' : 'O');
	89	Int_t nSec = (ring == 0 ? 20 : 40);
	90
	91	TObjArray* vtxArray = new TObjArray();
	92	vtxArray->SetName(Form("FMD%d%c",det,ringChar));
	93	detArray->AddAtAndExpand(vtxArray,ring);
	94	for(Int_t i = 0; i< nVtxbins; i++) {
	95	TH2F* hBg = pars->GetBackgroundCorrection(det, ringChar, i);
	96
	97	hMult = new TH2F(Form("FMD%d%c_vtxbin%d",det,ringChar,i),Form("FMD%d%c_vtxbin%d",det,ringChar,i),
	98	hBg->GetNbinsX(),
	99	hBg->GetXaxis()->GetXmin(),
	100	hBg->GetXaxis()->GetXmax(),
	101	nSec, 0, 2*TMath::Pi());
3bb122c7	102
bb8a464f	103	vtxArray->AddAtAndExpand(hMult,i);
3bb122c7	104	}
	105	}
	106	}
8dc7c4c2	107
bb8a464f	108
7c3e5162	109
3bb122c7	110
	111	}
	112	//_____________________________________________________________________
	113	void AliFMDAnalysisTaskDensity::ConnectInputData(Option_t /option*/)
	114	{
7c3e5162	115	if(fStandalone) {
	116	fESD = (AliESDFMD*)GetInputData(0);
	117	fVertex = (AliESDVertex*)GetInputData(1);
	118	}
3bb122c7	119	}
	120	//_____________________________________________________________________
	121	void AliFMDAnalysisTaskDensity::Exec(Option_t /option*/)
	122	{
	123	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
	124	AliFMDGeometry* geo = AliFMDGeometry::Instance();
	125
7c3e5162	126	//AliESDFMD* fmd = fESD->GetFMDData();
3bb122c7	127
3bb122c7	128	Double_t vertex[3];
7c3e5162	129	fVertex->GetXYZ(vertex);
8dc823cc	130	// Z Vtx cut
bb8a464f	131	if( TMath::Abs(vertex[2]) > pars->GetVtxCutZ()) {
bb8a464f	132	fStatus = kFALSE;
3bb122c7	133	return;
bb8a464f	134	}
	135	else
	136	fStatus = kTRUE;
	137
3bb122c7	138	Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
	139	Double_t vertexBinDouble = (vertex[2] + pars->GetVtxCutZ()) / delta;
	140
	141	Int_t vtxbin = (Int_t)vertexBinDouble;
1282ce49	142
8dc7c4c2	143	fVertexString.SetString(Form("%d",vtxbin));
8dc823cc	144	//Reset everything
3bb122c7	145	for(UShort_t det=1;det<=3;det++) {
8dc7c4c2	146	TObjArray* detArray = (TObjArray*)fArray.At(det);
3bb122c7	147	Int_t nRings = (det==1 ? 1 : 2);
	148	for (UShort_t ir = 0; ir < nRings; ir++) {
	149	TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
	150
	151	TH2F* hMult = (TH2F*)vtxArray->At(vtxbin);
	152	hMult->Reset();
	153	}
	154
	155	}
	156
	157
	158	for(UShort_t det=1;det<=3;det++) {
8dc7c4c2	159	TObjArray* detArray = (TObjArray*)fArray.At(det);
3bb122c7	160	Int_t nRings = (det==1 ? 1 : 2);
	161	for (UShort_t ir = 0; ir < nRings; ir++) {
	162	TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
	163
	164	TH2F* hMult = (TH2F*)vtxArray->At(vtxbin);
bb8a464f	165
3bb122c7	166	Char_t ring = (ir == 0 ? 'I' : 'O');
	167	UShort_t nsec = (ir == 0 ? 20 : 40);
	168	UShort_t nstr = (ir == 0 ? 512 : 256);
bb8a464f	169
3bb122c7	170	for(UShort_t sec =0; sec < nsec; sec++) {
3bb122c7	171	for(UShort_t strip = 0; strip < nstr; strip++) {
7c3e5162	172	Float_t mult = fESD->Multiplicity(det,ring,sec,strip);
6289b3e8	173	Float_t eta = fESD->Eta(det,ring,sec,strip);
6289b3e8	174	Float_t mult_cut = 0.2;
bb8a464f	175	if(mult == 0 \|\| mult == AliESDFMD::kInvalidMult) continue;
	176	//Particle number cut goes here...
	177	Float_t nParticles = 0;
6289b3e8	178	if(fESD->GetUniqueID() == kTRUE) {
	179	//proton + proton
	180	if(mult > mult_cut)
	181	nParticles = 1;
	182	}
	183	else {
	184
	185	//Pb+Pb
	186	Float_t mpv = pars->GetMPV(det,ring,eta);
	187	Float_t sigma = pars->GetSigma(det,ring,eta);
	188	Float_t alpha = pars->Get2MIPWeight(det,ring,eta);
	189	Float_t beta = pars->Get3MIPWeight(det,ring,eta);
	190
	191	Float_t sumCor = TMath::Landau(mult,mpv,sigma,kTRUE)+
	192	alphaTMath::Landau(mult,2mpv+2sigmaTMath::Log(2),2*sigma,kTRUE)+
	193	betaTMath::Landau(mult,3mpv+3sigmaTMath::Log(3),3*sigma,kTRUE);
	194	Float_t weight = TMath::Landau(mult,mpv,sigma,kTRUE)+
	195	2alphaTMath::Landau(mult,2mpv+2sigmaTMath::Log(2),2sigma,kTRUE)+
	196	3betaTMath::Landau(mult,3mpv+3sigmaTMath::Log(3),3sigma,kTRUE);
	197
	198
	199	if(mult > 0){//mult_cut) {
	200	if(sumCor) nParticles = weight / sumCor;
	201	else nParticles = 1;
	202	}
	203	//std::cout<<sumCor<<" "<<weight<<" "<<" "<<mult<<" "<<nParticles<<std::endl;
	204
	205	}
	206
	207
3bb122c7	208	Double_t x,y,z;
	209	geo->Detector2XYZ(det,ring,sec,strip,x,y,z);
	210	Float_t phi = TMath::ATan2(y,x);
	211	if(phi<0)
	212	phi = phi+2*TMath::Pi();
1282ce49	213	Float_t correction = GetAcceptanceCorrection(ring,strip);
6289b3e8	214	if(correction) nParticles = nParticles / correction;
bb8a464f	215	hMult->Fill(eta,phi,nParticles);
bb8a464f	216
3bb122c7	217
	218	}
	219	}
bb8a464f	220
3bb122c7	221	}
bb8a464f	222
3bb122c7	223
	224
	225	}
7c3e5162	226	if(fStandalone) {
	227	PostData(0, fOutputList);
	228	}
3bb122c7	229
	230	}
	231	//_____________________________________________________________________
1282ce49	232	Float_t AliFMDAnalysisTaskDensity::GetAcceptanceCorrection(Char_t ring, UShort_t strip)
	233	{
	234	AliFMDRing fmdring(ring);
	235	fmdring.Init();
	236	Float_t rad = fmdring.GetMaxR()-fmdring.GetMinR();
	237	Float_t segment = rad / fmdring.GetNStrips();
	238	Float_t radius = fmdring.GetMinR() + segment*strip;
	239
	240	Float_t basearea1 = 0.5fmdring.GetBaseStripLength(strip)TMath::Power(radius,2);
	241	Float_t basearea2 = 0.5fmdring.GetBaseStripLength(strip)TMath::Power((radius-segment),2);
	242	Float_t basearea = basearea1 - basearea2;
	243
	244	Float_t area1 = 0.5fmdring.GetStripLength(strip)TMath::Power(radius,2);
	245	Float_t area2 = 0.5fmdring.GetStripLength(strip)TMath::Power((radius-segment),2);
	246	Float_t area = area1 - area2;
	247
	248	Float_t correction = area/basearea;
	249
	250	return correction;
	251	}
	252
3bb122c7	253	//
	254	//EOF
	255	//