[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 "TF1.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 "AliFMDParameters.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());
  }
  
  fFuncPos = new TF1("funcPos","pol1",0,6);
  fFuncPos->SetParameters(0.99925,0.00298301);
  fFuncNeg = new TF1("funcNeg","pol1",-6,0);
  fFuncNeg->SetParameters(0.987583,-0.0101022);
  
  
}
//_____________________________________________________________________
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);
	  
	  if(mult == 0 || mult == AliESDFMD::kInvalidMult) continue;
	  //Particle number cut goes here...
	  //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 phi = pars->GetPhiFromSector(det,ring,sec);
	  Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]);
	  //std::cout<<phi<<"     "<<phicalc<<std::endl;
	  //  Float_t   r     = TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2));
	 // Float_t   theta = TMath::ATan2(r,z-vertex[2]);
	  // Float_t   etacalc   = -1*TMath::Log(TMath::Tan(0.5*theta));
	   
	   //  std::cout<<eta<<"    "<<etacalc<<std::endl;
	   //eta = etacalc;
	     
	  Float_t m   = pars->GetMPV(det,ring,eta);
	  Float_t s   = pars->GetSigma(det,ring,eta);
	  //AliFMDParameters* recopars = AliFMDParameters::Instance();
	  
	  Float_t mult_cut = 0.15;//m-2*s;//0.15;//0.2;//m-3*s;// 0.2;//0.01;//m-2*s;//0.2;
	  if(ring == 'I')
	    mult_cut = 0.10;
	  //mult_cut = (5*recopars->GetPedestalWidth(det,ring,sec,strip))/(recopars->GetPulseGain(det,ring,sec,strip)*recopars->GetDACPerMIP());
	  
	  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 > mult_cut) {
	      if(sumCor) nParticles = weight / sumCor;
	      else nParticles = 1;
	      
	    }
	    //std::cout<<sumCor<<"    "<<weight<<"    "<<"    "<<mult<<"  "<<nParticles<<std::endl;
	    
	  }
	  
	  
	  Float_t correction = GetAcceptanceCorrection(ring,strip);
	  
	  //std::cout<<"before "<<correction<<std::endl;
	  if(fESD->GetUniqueID() == kTRUE) {
	    if(det == 3) 
	      correction = correction / fFuncNeg->Eval(eta);
	    else
	      correction = correction / fFuncPos->Eval(eta);
	  }
	  
	  // std::cout<<correction<<std::endl;
	  if(correction) nParticles = nParticles / correction;
	  if(nParticles > 0)
	    hMult->Fill(eta,phi,nParticles);
	  
	  //if(det == 1 && ring =='I' && nParticles >0)
	  //if(nParticles > 0)
	  //  std::cout<<det<<"    "<<ring<<"    "<<sec<<"    "<<strip<<"   "<<mult<<std::endl;
	  
	}
      }
      
    }
    
	
  }
  

  if(fStandalone) {
    PostData(0, fOutputList); 
  }
  
}
//_____________________________________________________________________
Float_t AliFMDAnalysisTaskDensity::GetAcceptanceCorrection(Char_t ring, UShort_t strip)
{
  AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
  
  //AliFMDRing fmdring(ring);
  //fmdring.Init();
  Float_t   rad       = pars->GetMaxR(ring)-pars->GetMinR(ring);
  Float_t   nstrips   = (ring == 'I' ? 512 : 256);
  Float_t   segment   = rad / nstrips;
  Float_t   radius    = pars->GetMinR(ring) + segment*strip;
  
  Float_t   basearea1 = 0.5*pars->GetBaseStripLength(ring,strip)*TMath::Power(radius,2);
  Float_t   basearea2 = 0.5*pars->GetBaseStripLength(ring,strip)*TMath::Power((radius-segment),2);
  Float_t   basearea  = basearea1 - basearea2;
  
  Float_t   area1     = 0.5*pars->GetStripLength(ring,strip)*TMath::Power(radius,2);
  Float_t   area2     = 0.5*pars->GetStripLength(ring,strip)*TMath::Power((radius-segment),2);
  Float_t   area      = area1 - area2;
  
  Float_t correction = area/basearea;
  
  return correction;
}
//_____________________________________________________________________
/*Float_t AliFMDAnalysisTaskDensity::GetPhiFromSector(UShort_t det, Char_t ring, UShort_t sec)
{
  Int_t nsec = (ring == 'I' ? 20 : 40);
  Float_t basephi = 0;
  if(det == 1) 
    basephi = 1.72787594; 
  if(det == 2 && ring == 'I')
    basephi = 0.15707963;
  if(det == 2 && ring == 'O')
    basephi = 0.078539818;
  if(det == 3 && ring == 'I')
    basephi = 2.984513044;
  if(det == 3 && ring == 'O')
    basephi = 3.06305289;
  
  Float_t step = 2*TMath::Pi() / nsec;
  Float_t phi = 0;
  if(det == 3)
    phi = basephi - sec*step;
  else
    phi = basephi + sec*step;
  
  if(phi < 0) 
    phi = phi +2*TMath::Pi();
  if(phi > 2*TMath::Pi() )
    phi = phi - 2*TMath::Pi();
  
  return phi;
}

*/
//
//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"
f58a4769	11	#include "TF1.h"
3bb122c7	12	#include "AliFMDAnalysisTaskDensity.h"
	13	#include "AliAnalysisManager.h"
	14	#include "AliESDFMD.h"
	15	#include "AliESDEvent.h"
	16	#include "AliAODEvent.h"
	17	#include "AliAODHandler.h"
	18	#include "AliMCEventHandler.h"
	19	#include "AliStack.h"
	20	#include "AliESDVertex.h"
	21	#include "TMath.h"
	22	#include "AliFMDAnaParameters.h"
78f6f750	23	//#include "AliFMDParameters.h"
	24	//#include "AliFMDGeometry.h"
	25	//#include "AliFMDRing.h"
3bb122c7	26
	27	ClassImp(AliFMDAnalysisTaskDensity)
	28
	29	//_____________________________________________________________________
	30	AliFMDAnalysisTaskDensity::AliFMDAnalysisTaskDensity()
	31	: fDebug(0),
8dc7c4c2	32	fOutputList(),
8dc7c4c2	33	fArray(),
c78bc12b	34	fESD(0x0),
7c3e5162	35	fVertexString(),
7c3e5162	36	fVertex(0),
bb8a464f	37	fStandalone(kTRUE),
bb8a464f	38	fStatus(kTRUE)
3bb122c7	39	{
3bb122c7	40	// Default constructor
7c3e5162	41	DefineInput (0, AliESDFMD::Class());
7c3e5162	42	DefineInput (1, AliESDVertex::Class());
3bb122c7	43	DefineOutput(0,TList::Class());
	44	}
	45	//_____________________________________________________________________
7c3e5162	46	AliFMDAnalysisTaskDensity::AliFMDAnalysisTaskDensity(const char* name, Bool_t SE):
3bb122c7	47	AliAnalysisTask(name, "Density"),
3bb122c7	48	fDebug(0),
7c3e5162	49	fOutputList(0),
8dc7c4c2	50	fArray(),
c78bc12b	51	fESD(0x0),
7c3e5162	52	fVertexString(),
7c3e5162	53	fVertex(0),
bb8a464f	54	fStandalone(kTRUE),
bb8a464f	55	fStatus(kTRUE)
3bb122c7	56	{
7c3e5162	57	fStandalone = SE;
	58	if(fStandalone) {
	59	DefineInput (0, AliESDFMD::Class());
	60	DefineInput (1, AliESDVertex::Class());
	61	DefineOutput(0, TList::Class());
	62	}
f58a4769	63
	64	fFuncPos = new TF1("funcPos","pol1",0,6);
	65	fFuncPos->SetParameters(0.99925,0.00298301);
	66	fFuncNeg = new TF1("funcNeg","pol1",-6,0);
	67	fFuncNeg->SetParameters(0.987583,-0.0101022);
78f6f750	68
f58a4769	69
3bb122c7	70	}
	71	//_____________________________________________________________________
	72	void AliFMDAnalysisTaskDensity::CreateOutputObjects()
	73	{
	74	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
3bb122c7	75
8dc7c4c2	76	fArray.SetName("FMD");
8dc7c4c2	77	fArray.SetOwner();
7c3e5162	78	if(!fOutputList)
	79	fOutputList = new TList();
	80	fOutputList->SetName("density_list");
3bb122c7	81
bb8a464f	82	fOutputList->Add(&fArray);
	83	fOutputList->Add(&fVertexString);
	84
3bb122c7	85	TH2F* hMult = 0;
	86
	87	Int_t nVtxbins = pars->GetNvtxBins();
	88
	89	for(Int_t det =1; det<=3;det++)
	90	{
	91	TObjArray* detArray = new TObjArray();
	92	detArray->SetName(Form("FMD%d",det));
8dc7c4c2	93	fArray.AddAtAndExpand(detArray,det);
3bb122c7	94	Int_t nRings = (det==1 ? 1 : 2);
	95	for(Int_t ring = 0;ring<nRings;ring++)
	96	{
	97	Char_t ringChar = (ring == 0 ? 'I' : 'O');
	98	Int_t nSec = (ring == 0 ? 20 : 40);
	99
	100	TObjArray* vtxArray = new TObjArray();
	101	vtxArray->SetName(Form("FMD%d%c",det,ringChar));
	102	detArray->AddAtAndExpand(vtxArray,ring);
	103	for(Int_t i = 0; i< nVtxbins; i++) {
	104	TH2F* hBg = pars->GetBackgroundCorrection(det, ringChar, i);
	105
	106	hMult = new TH2F(Form("FMD%d%c_vtxbin%d",det,ringChar,i),Form("FMD%d%c_vtxbin%d",det,ringChar,i),
	107	hBg->GetNbinsX(),
	108	hBg->GetXaxis()->GetXmin(),
	109	hBg->GetXaxis()->GetXmax(),
	110	nSec, 0, 2*TMath::Pi());
3bb122c7	111
bb8a464f	112	vtxArray->AddAtAndExpand(hMult,i);
3bb122c7	113	}
	114	}
	115	}
8dc7c4c2	116
bb8a464f	117
7c3e5162	118
3bb122c7	119
	120	}
	121	//_____________________________________________________________________
	122	void AliFMDAnalysisTaskDensity::ConnectInputData(Option_t /option*/)
	123	{
7c3e5162	124	if(fStandalone) {
	125	fESD = (AliESDFMD*)GetInputData(0);
	126	fVertex = (AliESDVertex*)GetInputData(1);
	127	}
3bb122c7	128	}
f58a4769	129
	130
	131
	132
3bb122c7	133	//_____________________________________________________________________
	134	void AliFMDAnalysisTaskDensity::Exec(Option_t /option*/)
	135	{
	136	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
78f6f750	137	// AliFMDGeometry* geo = AliFMDGeometry::Instance();
3bb122c7	138
7c3e5162	139	//AliESDFMD* fmd = fESD->GetFMDData();
3bb122c7	140
3bb122c7	141	Double_t vertex[3];
7c3e5162	142	fVertex->GetXYZ(vertex);
8dc823cc	143	// Z Vtx cut
bb8a464f	144	if( TMath::Abs(vertex[2]) > pars->GetVtxCutZ()) {
bb8a464f	145	fStatus = kFALSE;
3bb122c7	146	return;
bb8a464f	147	}
	148	else
	149	fStatus = kTRUE;
	150
3bb122c7	151	Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
	152	Double_t vertexBinDouble = (vertex[2] + pars->GetVtxCutZ()) / delta;
	153
	154	Int_t vtxbin = (Int_t)vertexBinDouble;
f58a4769	155
1282ce49	156
f58a4769	157
8dc7c4c2	158	fVertexString.SetString(Form("%d",vtxbin));
8dc823cc	159	//Reset everything
3bb122c7	160	for(UShort_t det=1;det<=3;det++) {
8dc7c4c2	161	TObjArray* detArray = (TObjArray*)fArray.At(det);
3bb122c7	162	Int_t nRings = (det==1 ? 1 : 2);
	163	for (UShort_t ir = 0; ir < nRings; ir++) {
	164	TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
	165
	166	TH2F* hMult = (TH2F*)vtxArray->At(vtxbin);
	167	hMult->Reset();
	168	}
	169
	170	}
	171
	172
	173	for(UShort_t det=1;det<=3;det++) {
8dc7c4c2	174	TObjArray* detArray = (TObjArray*)fArray.At(det);
3bb122c7	175	Int_t nRings = (det==1 ? 1 : 2);
	176	for (UShort_t ir = 0; ir < nRings; ir++) {
	177	TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
	178
	179	TH2F* hMult = (TH2F*)vtxArray->At(vtxbin);
bb8a464f	180
3bb122c7	181	Char_t ring = (ir == 0 ? 'I' : 'O');
	182	UShort_t nsec = (ir == 0 ? 20 : 40);
	183	UShort_t nstr = (ir == 0 ? 512 : 256);
bb8a464f	184
3bb122c7	185	for(UShort_t sec =0; sec < nsec; sec++) {
3bb122c7	186	for(UShort_t strip = 0; strip < nstr; strip++) {
7c3e5162	187	Float_t mult = fESD->Multiplicity(det,ring,sec,strip);
f58a4769	188	//Float_t eta = fESD->Eta(det,ring,sec,strip);
cc066cb9	189
bb8a464f	190	if(mult == 0 \|\| mult == AliESDFMD::kInvalidMult) continue;
bb8a464f	191	//Particle number cut goes here...
f58a4769	192	//Double_t x,y,z;
	193	//geo->Detector2XYZ(det,ring,sec,strip,x,y,z);
	194	// Float_t phi = TMath::ATan2(y,x);
	195	// if(phi<0)
	196	// phi = phi+2*TMath::Pi();
	197
	198	Float_t phi = pars->GetPhiFromSector(det,ring,sec);
	199	Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]);
	200	//std::cout<<phi<<" "<<phicalc<<std::endl;
	201	// Float_t r = TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2));
	202	// Float_t theta = TMath::ATan2(r,z-vertex[2]);
	203	// Float_t etacalc = -1TMath::Log(TMath::Tan(0.5theta));
	204
	205	// std::cout<<eta<<" "<<etacalc<<std::endl;
	206	//eta = etacalc;
	207
cc066cb9	208	Float_t m = pars->GetMPV(det,ring,eta);
cc066cb9	209	Float_t s = pars->GetSigma(det,ring,eta);
f58a4769	210	//AliFMDParameters* recopars = AliFMDParameters::Instance();
cc066cb9	211
f58a4769	212	Float_t mult_cut = 0.15;//m-2s;//0.15;//0.2;//m-3s;// 0.2;//0.01;//m-2*s;//0.2;
	213	if(ring == 'I')
	214	mult_cut = 0.10;
	215	//mult_cut = (5recopars->GetPedestalWidth(det,ring,sec,strip))/(recopars->GetPulseGain(det,ring,sec,strip)recopars->GetDACPerMIP());
cc066cb9	216
bb8a464f	217	Float_t nParticles = 0;
6289b3e8	218	if(fESD->GetUniqueID() == kTRUE) {
6289b3e8	219	//proton + proton
cc066cb9	220
cc066cb9	221	if(mult > mult_cut) {
6289b3e8	222	nParticles = 1;
cc066cb9	223
cc066cb9	224	}
6289b3e8	225	}
	226	else {
	227
	228	//Pb+Pb
	229	Float_t mpv = pars->GetMPV(det,ring,eta);
	230	Float_t sigma = pars->GetSigma(det,ring,eta);
	231	Float_t alpha = pars->Get2MIPWeight(det,ring,eta);
	232	Float_t beta = pars->Get3MIPWeight(det,ring,eta);
	233
	234	Float_t sumCor = TMath::Landau(mult,mpv,sigma,kTRUE)+
	235	alphaTMath::Landau(mult,2mpv+2sigmaTMath::Log(2),2*sigma,kTRUE)+
	236	betaTMath::Landau(mult,3mpv+3sigmaTMath::Log(3),3*sigma,kTRUE);
	237	Float_t weight = TMath::Landau(mult,mpv,sigma,kTRUE)+
	238	2alphaTMath::Landau(mult,2mpv+2sigmaTMath::Log(2),2sigma,kTRUE)+
	239	3betaTMath::Landau(mult,3mpv+3sigmaTMath::Log(3),3sigma,kTRUE);
	240
	241
cc066cb9	242	if(mult > mult_cut) {
6289b3e8	243	if(sumCor) nParticles = weight / sumCor;
6289b3e8	244	else nParticles = 1;
cc066cb9	245
6289b3e8	246	}
	247	//std::cout<<sumCor<<" "<<weight<<" "<<" "<<mult<<" "<<nParticles<<std::endl;
	248
	249	}
	250
	251
cc066cb9	252
	253
	254
1282ce49	255	Float_t correction = GetAcceptanceCorrection(ring,strip);
f58a4769	256
f58a4769	257	//std::cout<<"before "<<correction<<std::endl;
78f6f750	258	if(fESD->GetUniqueID() == kTRUE) {
	259	if(det == 3)
	260	correction = correction / fFuncNeg->Eval(eta);
	261	else
	262	correction = correction / fFuncPos->Eval(eta);
	263	}
f58a4769	264
f58a4769	265	// std::cout<<correction<<std::endl;
6289b3e8	266	if(correction) nParticles = nParticles / correction;
cc066cb9	267	if(nParticles > 0)
cc066cb9	268	hMult->Fill(eta,phi,nParticles);
bb8a464f	269
cc066cb9	270	//if(det == 1 && ring =='I' && nParticles >0)
f58a4769	271	//if(nParticles > 0)
f58a4769	272	// std::cout<<det<<" "<<ring<<" "<<sec<<" "<<strip<<" "<<mult<<std::endl;
3bb122c7	273
	274	}
	275	}
bb8a464f	276
3bb122c7	277	}
bb8a464f	278
3bb122c7	279
	280
	281	}
78f6f750	282
78f6f750	283
7c3e5162	284	if(fStandalone) {
	285	PostData(0, fOutputList);
	286	}
3bb122c7	287
	288	}
	289	//_____________________________________________________________________
1282ce49	290	Float_t AliFMDAnalysisTaskDensity::GetAcceptanceCorrection(Char_t ring, UShort_t strip)
1282ce49	291	{
78f6f750	292	AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
1282ce49	293
78f6f750	294	//AliFMDRing fmdring(ring);
	295	//fmdring.Init();
	296	Float_t rad = pars->GetMaxR(ring)-pars->GetMinR(ring);
	297	Float_t nstrips = (ring == 'I' ? 512 : 256);
	298	Float_t segment = rad / nstrips;
	299	Float_t radius = pars->GetMinR(ring) + segment*strip;
	300
	301	Float_t basearea1 = 0.5pars->GetBaseStripLength(ring,strip)TMath::Power(radius,2);
	302	Float_t basearea2 = 0.5pars->GetBaseStripLength(ring,strip)TMath::Power((radius-segment),2);
1282ce49	303	Float_t basearea = basearea1 - basearea2;
1282ce49	304
78f6f750	305	Float_t area1 = 0.5pars->GetStripLength(ring,strip)TMath::Power(radius,2);
78f6f750	306	Float_t area2 = 0.5pars->GetStripLength(ring,strip)TMath::Power((radius-segment),2);
1282ce49	307	Float_t area = area1 - area2;
	308
	309	Float_t correction = area/basearea;
	310
	311	return correction;
	312	}
f58a4769	313	//_____________________________________________________________________
78f6f750	314	/*Float_t AliFMDAnalysisTaskDensity::GetPhiFromSector(UShort_t det, Char_t ring, UShort_t sec)
f58a4769	315	{
	316	Int_t nsec = (ring == 'I' ? 20 : 40);
	317	Float_t basephi = 0;
	318	if(det == 1)
	319	basephi = 1.72787594;
	320	if(det == 2 && ring == 'I')
	321	basephi = 0.15707963;
	322	if(det == 2 && ring == 'O')
	323	basephi = 0.078539818;
	324	if(det == 3 && ring == 'I')
	325	basephi = 2.984513044;
	326	if(det == 3 && ring == 'O')
	327	basephi = 3.06305289;
	328
	329	Float_t step = 2*TMath::Pi() / nsec;
	330	Float_t phi = 0;
	331	if(det == 3)
	332	phi = basephi - sec*step;
	333	else
	334	phi = basephi + sec*step;
	335
	336	if(phi < 0)
	337	phi = phi +2*TMath::Pi();
	338	if(phi > 2*TMath::Pi() )
	339	phi = phi - 2*TMath::Pi();
	340
	341	return phi;
	342	}
	343
78f6f750	344	*/
3bb122c7	345	//
	346	//EOF
	347	//