]>
Commit | Line | Data |
---|---|---|
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 "TF1.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" | |
19 | #include "AliStack.h" | |
20 | #include "AliESDVertex.h" | |
21 | #include "TMath.h" | |
22 | #include "AliFMDAnaParameters.h" | |
23 | //#include "AliFMDParameters.h" | |
24 | //#include "AliFMDGeometry.h" | |
25 | //#include "AliFMDRing.h" | |
26 | ||
27 | ClassImp(AliFMDAnalysisTaskDensity) | |
28 | ||
29 | //_____________________________________________________________________ | |
30 | AliFMDAnalysisTaskDensity::AliFMDAnalysisTaskDensity() | |
31 | : fDebug(0), | |
32 | fOutputList(), | |
33 | fESD(0x0), | |
34 | fVertexString(), | |
35 | fVertex(0), | |
36 | fStandalone(kTRUE), | |
37 | fStatus(kTRUE) | |
38 | { | |
39 | // Default constructor | |
40 | DefineInput (0, AliESDFMD::Class()); | |
41 | DefineInput (1, AliESDVertex::Class()); | |
42 | DefineOutput(0,TList::Class()); | |
43 | } | |
44 | //_____________________________________________________________________ | |
45 | AliFMDAnalysisTaskDensity::AliFMDAnalysisTaskDensity(const char* name, Bool_t SE): | |
46 | AliAnalysisTask(name, "Density"), | |
47 | fDebug(0), | |
48 | fOutputList(0), | |
49 | fESD(0x0), | |
50 | fVertexString(), | |
51 | fVertex(0), | |
52 | fStandalone(kTRUE), | |
53 | fStatus(kTRUE) | |
54 | { | |
55 | fStandalone = SE; | |
56 | if(fStandalone) { | |
57 | DefineInput (0, AliESDFMD::Class()); | |
58 | DefineInput (1, AliESDVertex::Class()); | |
59 | DefineOutput(0, TList::Class()); | |
60 | } | |
61 | ||
62 | } | |
63 | //_____________________________________________________________________ | |
64 | void AliFMDAnalysisTaskDensity::CreateOutputObjects() | |
65 | { | |
66 | AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance(); | |
67 | ||
68 | if(!fOutputList) | |
69 | fOutputList = new TList(); | |
70 | fOutputList->SetName("density_list"); | |
71 | ||
72 | fOutputList->Add(&fVertexString); | |
73 | ||
74 | TH2F* hMult = 0; | |
75 | ||
76 | Int_t nVtxbins = pars->GetNvtxBins(); | |
77 | ||
78 | for(Int_t det =1; det<=3;det++) | |
79 | { | |
80 | Int_t nRings = (det==1 ? 1 : 2); | |
81 | for(Int_t ring = 0;ring<nRings;ring++) | |
82 | { | |
83 | Char_t ringChar = (ring == 0 ? 'I' : 'O'); | |
84 | Int_t nSec = (ring == 0 ? 20 : 40); | |
85 | ||
86 | for(Int_t i = 0; i< nVtxbins; i++) { | |
87 | TH2F* hBg = pars->GetBackgroundCorrection(det, ringChar, i); | |
88 | ||
89 | hMult = new TH2F(Form("FMD%d%c_vtxbin%d",det,ringChar,i),Form("FMD%d%c_vtxbin%d",det,ringChar,i), | |
90 | hBg->GetNbinsX(), | |
91 | hBg->GetXaxis()->GetXmin(), | |
92 | hBg->GetXaxis()->GetXmax(), | |
93 | nSec, 0, 2*TMath::Pi()); | |
94 | ||
95 | fOutputList->Add(hMult); | |
96 | } | |
97 | } | |
98 | } | |
99 | ||
100 | ||
101 | ||
102 | ||
103 | } | |
104 | //_____________________________________________________________________ | |
105 | void AliFMDAnalysisTaskDensity::ConnectInputData(Option_t */*option*/) | |
106 | { | |
107 | if(fStandalone) { | |
108 | fESD = (AliESDFMD*)GetInputData(0); | |
109 | fVertex = (AliESDVertex*)GetInputData(1); | |
110 | } | |
111 | } | |
112 | ||
113 | ||
114 | ||
115 | ||
116 | //_____________________________________________________________________ | |
117 | void AliFMDAnalysisTaskDensity::Exec(Option_t */*option*/) | |
118 | { | |
119 | AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance(); | |
120 | // AliFMDGeometry* geo = AliFMDGeometry::Instance(); | |
121 | ||
122 | //AliESDFMD* fmd = fESD->GetFMDData(); | |
123 | ||
124 | Double_t vertex[3]; | |
125 | fVertex->GetXYZ(vertex); | |
126 | // Z Vtx cut | |
127 | if( TMath::Abs(vertex[2]) > pars->GetVtxCutZ()) { | |
128 | fStatus = kFALSE; | |
129 | return; | |
130 | } | |
131 | else | |
132 | fStatus = kTRUE; | |
133 | ||
134 | Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins(); | |
135 | Double_t vertexBinDouble = (vertex[2] + pars->GetVtxCutZ()) / delta; | |
136 | ||
137 | Int_t vtxbin = (Int_t)vertexBinDouble; | |
138 | ||
139 | fVertexString.SetString(Form("%d",vtxbin)); | |
140 | //Reset everything | |
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)); | |
146 | hMult->Reset(); | |
147 | } | |
148 | ||
149 | } | |
150 | ||
151 | ||
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++) { | |
155 | ||
156 | Char_t ring = (ir == 0 ? 'I' : 'O'); | |
157 | TH2F* hMult = (TH2F*)fOutputList->FindObject(Form("FMD%d%c_vtxbin%d",det,ring,vtxbin)); | |
158 | ||
159 | UShort_t nsec = (ir == 0 ? 20 : 40); | |
160 | UShort_t nstr = (ir == 0 ? 512 : 256); | |
161 | ||
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); | |
165 | ||
166 | if(mult == 0 || mult == AliESDFMD::kInvalidMult) continue; | |
167 | ||
168 | Float_t phi = pars->GetPhiFromSector(det,ring,sec); | |
169 | Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]); | |
170 | ||
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; | |
172 | if(ring == 'I') | |
173 | mult_cut = 0.10; | |
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) { | |
177 | //proton + proton | |
178 | ||
179 | if(mult > mult_cut) { | |
180 | nParticles = 1; | |
181 | ||
182 | } | |
183 | } | |
184 | else { | |
185 | ||
186 | //Pb+Pb | |
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); | |
191 | ||
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); | |
198 | ||
199 | ||
200 | if(mult > mult_cut) { | |
201 | if(sumCor) nParticles = weight / sumCor; | |
202 | else nParticles = 1; | |
203 | ||
204 | } | |
205 | //std::cout<<sumCor<<" "<<weight<<" "<<" "<<mult<<" "<<nParticles<<std::endl; | |
206 | ||
207 | } | |
208 | ||
209 | ||
210 | ||
211 | ||
212 | ||
213 | Float_t correction = GetAcceptanceCorrection(ring,strip); | |
214 | ||
215 | //std::cout<<"before "<<correction<<std::endl; | |
216 | if(fESD->GetUniqueID() == kTRUE) { | |
217 | TH1F* hDoubleHitCorrection = pars->GetDoubleHitCorrection(det,ring); | |
218 | ||
219 | if(hDoubleHitCorrection->GetBinContent(hDoubleHitCorrection->FindBin(eta)) != 0) | |
220 | correction = correction*hDoubleHitCorrection->GetBinContent(hDoubleHitCorrection->FindBin(eta)); | |
221 | ||
222 | } | |
223 | ||
224 | if(correction) nParticles = nParticles / correction; | |
225 | if(nParticles > 0) | |
226 | hMult->Fill(eta,phi,nParticles); | |
227 | ||
228 | ||
229 | } | |
230 | } | |
231 | ||
232 | } | |
233 | ||
234 | ||
235 | ||
236 | } | |
237 | ||
238 | ||
239 | if(fStandalone) { | |
240 | PostData(0, fOutputList); | |
241 | } | |
242 | ||
243 | } | |
244 | //_____________________________________________________________________ | |
245 | Float_t AliFMDAnalysisTaskDensity::GetAcceptanceCorrection(Char_t ring, UShort_t strip) | |
246 | { | |
247 | AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance(); | |
248 | ||
249 | //AliFMDRing fmdring(ring); | |
250 | //fmdring.Init(); | |
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; | |
255 | ||
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; | |
259 | ||
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; | |
263 | ||
264 | Float_t correction = area/basearea; | |
265 | ||
266 | return correction; | |
267 | } | |
268 | //_____________________________________________________________________ | |
269 | // | |
270 | //EOF | |
271 | // |