4 #include <TInterpreter.h>
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"
20 #include "AliESDVertex.h"
22 #include "AliFMDAnaParameters.h"
23 //#include "AliFMDParameters.h"
24 //#include "AliFMDGeometry.h"
25 //#include "AliFMDRing.h"
27 ClassImp(AliFMDAnalysisTaskDensity)
29 //_____________________________________________________________________
30 AliFMDAnalysisTaskDensity::AliFMDAnalysisTaskDensity()
39 // Default constructor
40 DefineInput (0, AliESDFMD::Class());
41 DefineInput (1, AliESDVertex::Class());
42 DefineOutput(0,TList::Class());
44 //_____________________________________________________________________
45 AliFMDAnalysisTaskDensity::AliFMDAnalysisTaskDensity(const char* name, Bool_t SE):
46 AliAnalysisTask(name, "Density"),
57 DefineInput (0, AliESDFMD::Class());
58 DefineInput (1, AliESDVertex::Class());
59 DefineOutput(0, TList::Class());
63 //_____________________________________________________________________
64 void AliFMDAnalysisTaskDensity::CreateOutputObjects()
66 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
69 fOutputList = new TList();
70 fOutputList->SetName("density_list");
72 fOutputList->Add(&fVertexString);
76 Int_t nVtxbins = pars->GetNvtxBins();
78 for(Int_t det =1; det<=3;det++)
80 Int_t nRings = (det==1 ? 1 : 2);
81 for(Int_t ring = 0;ring<nRings;ring++)
83 Char_t ringChar = (ring == 0 ? 'I' : 'O');
84 Int_t nSec = (ring == 0 ? 20 : 40);
86 for(Int_t i = 0; i< nVtxbins; i++) {
87 TH2F* hBg = pars->GetBackgroundCorrection(det, ringChar, i);
89 hMult = new TH2F(Form("FMD%d%c_vtxbin%d",det,ringChar,i),Form("FMD%d%c_vtxbin%d",det,ringChar,i),
91 hBg->GetXaxis()->GetXmin(),
92 hBg->GetXaxis()->GetXmax(),
93 nSec, 0, 2*TMath::Pi());
95 fOutputList->Add(hMult);
104 //_____________________________________________________________________
105 void AliFMDAnalysisTaskDensity::ConnectInputData(Option_t */*option*/)
108 fESD = (AliESDFMD*)GetInputData(0);
109 fVertex = (AliESDVertex*)GetInputData(1);
116 //_____________________________________________________________________
117 void AliFMDAnalysisTaskDensity::Exec(Option_t */*option*/)
119 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
120 // AliFMDGeometry* geo = AliFMDGeometry::Instance();
122 //AliESDFMD* fmd = fESD->GetFMDData();
125 fVertex->GetXYZ(vertex);
127 if( TMath::Abs(vertex[2]) > pars->GetVtxCutZ()) {
134 Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
135 Double_t vertexBinDouble = (vertex[2] + pars->GetVtxCutZ()) / delta;
137 Int_t vtxbin = (Int_t)vertexBinDouble;
139 fVertexString.SetString(Form("%d",vtxbin));
141 for(UShort_t det=1;det<=3;det++) {
142 Int_t nRings = (det==1 ? 1 : 2);
143 for (UShort_t ir = 0; ir < nRings; ir++) {
144 Char_t ring = (ir == 0 ? 'I' : 'O');
145 TH2F* hMult = (TH2F*)fOutputList->FindObject(Form("FMD%d%c_vtxbin%d",det,ring,vtxbin));
152 for(UShort_t det=1;det<=3;det++) {
153 Int_t nRings = (det==1 ? 1 : 2);
154 for (UShort_t ir = 0; ir < nRings; ir++) {
156 Char_t ring = (ir == 0 ? 'I' : 'O');
157 TH2F* hMult = (TH2F*)fOutputList->FindObject(Form("FMD%d%c_vtxbin%d",det,ring,vtxbin));
159 UShort_t nsec = (ir == 0 ? 20 : 40);
160 UShort_t nstr = (ir == 0 ? 512 : 256);
162 for(UShort_t sec =0; sec < nsec; sec++) {
163 for(UShort_t strip = 0; strip < nstr; strip++) {
164 Float_t mult = fESD->Multiplicity(det,ring,sec,strip);
166 if(mult == 0 || mult == AliESDFMD::kInvalidMult) continue;
168 Float_t phi = pars->GetPhiFromSector(det,ring,sec);
169 Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]);
171 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;
174 //Float_t mult_cut = pars->GetMPV(det,ring,eta) - 5*pars->GetSigma(det,ring,eta);
175 Float_t nParticles = 0;
176 if(fESD->GetUniqueID() == kTRUE) {
179 if(mult > mult_cut) {
187 Float_t mpv = pars->GetMPV(det,ring,eta);
188 Float_t sigma = pars->GetSigma(det,ring,eta);
189 Float_t alpha = pars->Get2MIPWeight(det,ring,eta);
190 Float_t beta = pars->Get3MIPWeight(det,ring,eta);
192 Float_t sumCor = TMath::Landau(mult,mpv,sigma,kTRUE)+
193 alpha*TMath::Landau(mult,2*mpv+2*sigma*TMath::Log(2),2*sigma,kTRUE)+
194 beta*TMath::Landau(mult,3*mpv+3*sigma*TMath::Log(3),3*sigma,kTRUE);
195 Float_t weight = TMath::Landau(mult,mpv,sigma,kTRUE)+
196 2*alpha*TMath::Landau(mult,2*mpv+2*sigma*TMath::Log(2),2*sigma,kTRUE)+
197 3*beta*TMath::Landau(mult,3*mpv+3*sigma*TMath::Log(3),3*sigma,kTRUE);
200 if(mult > mult_cut) {
201 if(sumCor) nParticles = weight / sumCor;
205 //std::cout<<sumCor<<" "<<weight<<" "<<" "<<mult<<" "<<nParticles<<std::endl;
213 Float_t correction = GetAcceptanceCorrection(ring,strip);
215 //std::cout<<"before "<<correction<<std::endl;
216 if(fESD->GetUniqueID() == kTRUE) {
217 TH1F* hDoubleHitCorrection = pars->GetDoubleHitCorrection(det,ring);
219 if(hDoubleHitCorrection->GetBinContent(hDoubleHitCorrection->FindBin(eta)) != 0)
220 correction = correction*hDoubleHitCorrection->GetBinContent(hDoubleHitCorrection->FindBin(eta));
224 if(correction) nParticles = nParticles / correction;
226 hMult->Fill(eta,phi,nParticles);
240 PostData(0, fOutputList);
244 //_____________________________________________________________________
245 Float_t AliFMDAnalysisTaskDensity::GetAcceptanceCorrection(Char_t ring, UShort_t strip)
247 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
249 //AliFMDRing fmdring(ring);
251 Float_t rad = pars->GetMaxR(ring)-pars->GetMinR(ring);
252 Float_t nstrips = (ring == 'I' ? 512 : 256);
253 Float_t segment = rad / nstrips;
254 Float_t radius = pars->GetMinR(ring) + segment*strip;
256 Float_t basearea1 = 0.5*pars->GetBaseStripLength(ring,strip)*TMath::Power(radius,2);
257 Float_t basearea2 = 0.5*pars->GetBaseStripLength(ring,strip)*TMath::Power((radius-segment),2);
258 Float_t basearea = basearea1 - basearea2;
260 Float_t area1 = 0.5*pars->GetStripLength(ring,strip)*TMath::Power(radius,2);
261 Float_t area2 = 0.5*pars->GetStripLength(ring,strip)*TMath::Power((radius-segment),2);
262 Float_t area = area1 - area2;
264 Float_t correction = area/basearea;
268 //_____________________________________________________________________