]>
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" | |
da0805e2 | 20 | #include "AliFMDAnaCalibEnergyDistribution.h" |
3bb122c7 | 21 | #include "AliESDVertex.h" |
22 | #include "TMath.h" | |
23 | #include "AliFMDAnaParameters.h" | |
78f6f750 | 24 | //#include "AliFMDParameters.h" |
25 | //#include "AliFMDGeometry.h" | |
26 | //#include "AliFMDRing.h" | |
3bb122c7 | 27 | |
28 | ClassImp(AliFMDAnalysisTaskDensity) | |
29 | ||
30 | //_____________________________________________________________________ | |
31 | AliFMDAnalysisTaskDensity::AliFMDAnalysisTaskDensity() | |
32 | : fDebug(0), | |
8dc7c4c2 | 33 | fOutputList(), |
c78bc12b | 34 | fESD(0x0), |
7c3e5162 | 35 | fVertexString(), |
36 | fVertex(0), | |
bb8a464f | 37 | fStandalone(kTRUE), |
38 | fStatus(kTRUE) | |
3bb122c7 | 39 | { |
40 | // Default constructor | |
7c3e5162 | 41 | DefineInput (0, AliESDFMD::Class()); |
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"), |
48 | fDebug(0), | |
7c3e5162 | 49 | fOutputList(0), |
c78bc12b | 50 | fESD(0x0), |
7c3e5162 | 51 | fVertexString(), |
52 | fVertex(0), | |
bb8a464f | 53 | fStandalone(kTRUE), |
54 | fStatus(kTRUE) | |
3bb122c7 | 55 | { |
7c3e5162 | 56 | fStandalone = SE; |
57 | if(fStandalone) { | |
58 | DefineInput (0, AliESDFMD::Class()); | |
59 | DefineInput (1, AliESDVertex::Class()); | |
60 | DefineOutput(0, TList::Class()); | |
61 | } | |
f58a4769 | 62 | |
3bb122c7 | 63 | } |
64 | //_____________________________________________________________________ | |
65 | void AliFMDAnalysisTaskDensity::CreateOutputObjects() | |
66 | { | |
67 | AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance(); | |
3bb122c7 | 68 | |
da0805e2 | 69 | |
70 | ||
7c3e5162 | 71 | if(!fOutputList) |
72 | fOutputList = new TList(); | |
73 | fOutputList->SetName("density_list"); | |
3bb122c7 | 74 | |
bb8a464f | 75 | fOutputList->Add(&fVertexString); |
76 | ||
f55d559b | 77 | |
78 | ||
3bb122c7 | 79 | TH2F* hMult = 0; |
80 | ||
81 | Int_t nVtxbins = pars->GetNvtxBins(); | |
82 | ||
f55d559b | 83 | |
3bb122c7 | 84 | for(Int_t det =1; det<=3;det++) |
85 | { | |
3bb122c7 | 86 | Int_t nRings = (det==1 ? 1 : 2); |
87 | for(Int_t ring = 0;ring<nRings;ring++) | |
88 | { | |
89 | Char_t ringChar = (ring == 0 ? 'I' : 'O'); | |
90 | Int_t nSec = (ring == 0 ? 20 : 40); | |
91 | ||
3bb122c7 | 92 | for(Int_t i = 0; i< nVtxbins; i++) { |
93 | TH2F* hBg = pars->GetBackgroundCorrection(det, ringChar, i); | |
94 | ||
95 | hMult = new TH2F(Form("FMD%d%c_vtxbin%d",det,ringChar,i),Form("FMD%d%c_vtxbin%d",det,ringChar,i), | |
96 | hBg->GetNbinsX(), | |
97 | hBg->GetXaxis()->GetXmin(), | |
98 | hBg->GetXaxis()->GetXmax(), | |
99 | nSec, 0, 2*TMath::Pi()); | |
f7356393 | 100 | hMult->Sumw2(); |
3bb122c7 | 101 | |
9f55be54 | 102 | fOutputList->Add(hMult); |
3bb122c7 | 103 | } |
104 | } | |
105 | } | |
8dc7c4c2 | 106 | |
bb8a464f | 107 | |
7c3e5162 | 108 | |
3bb122c7 | 109 | |
110 | } | |
111 | //_____________________________________________________________________ | |
112 | void AliFMDAnalysisTaskDensity::ConnectInputData(Option_t */*option*/) | |
113 | { | |
7c3e5162 | 114 | if(fStandalone) { |
115 | fESD = (AliESDFMD*)GetInputData(0); | |
116 | fVertex = (AliESDVertex*)GetInputData(1); | |
117 | } | |
3bb122c7 | 118 | } |
da0805e2 | 119 | //_____________________________________________________________________ |
120 | void AliFMDAnalysisTaskDensity::Init() { | |
121 | //TFile f("/home/canute/ALICE/FMDanalysis/FirstAnalysis/energydistributions_0_0_1_0_0_0.root"); | |
122 | //fEnergyDistribution = dynamic_cast<AliFMDAnaCalibEnergyDistribution*>(f.Get("energydistributions")); | |
f58a4769 | 123 | |
da0805e2 | 124 | } |
3bb122c7 | 125 | //_____________________________________________________________________ |
126 | void AliFMDAnalysisTaskDensity::Exec(Option_t */*option*/) | |
127 | { | |
128 | AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance(); | |
78f6f750 | 129 | // AliFMDGeometry* geo = AliFMDGeometry::Instance(); |
3bb122c7 | 130 | |
7c3e5162 | 131 | //AliESDFMD* fmd = fESD->GetFMDData(); |
3bb122c7 | 132 | |
da0805e2 | 133 | |
3bb122c7 | 134 | Double_t vertex[3]; |
7c3e5162 | 135 | fVertex->GetXYZ(vertex); |
8dc823cc | 136 | // Z Vtx cut |
da0805e2 | 137 | |
bb8a464f | 138 | if( TMath::Abs(vertex[2]) > pars->GetVtxCutZ()) { |
139 | fStatus = kFALSE; | |
3bb122c7 | 140 | return; |
bb8a464f | 141 | } |
142 | else | |
143 | fStatus = kTRUE; | |
144 | ||
3bb122c7 | 145 | Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins(); |
146 | Double_t vertexBinDouble = (vertex[2] + pars->GetVtxCutZ()) / delta; | |
147 | ||
148 | Int_t vtxbin = (Int_t)vertexBinDouble; | |
f58a4769 | 149 | |
8dc7c4c2 | 150 | fVertexString.SetString(Form("%d",vtxbin)); |
8dc823cc | 151 | //Reset everything |
3bb122c7 | 152 | for(UShort_t det=1;det<=3;det++) { |
3bb122c7 | 153 | Int_t nRings = (det==1 ? 1 : 2); |
154 | for (UShort_t ir = 0; ir < nRings; ir++) { | |
9f55be54 | 155 | Char_t ring = (ir == 0 ? 'I' : 'O'); |
156 | TH2F* hMult = (TH2F*)fOutputList->FindObject(Form("FMD%d%c_vtxbin%d",det,ring,vtxbin)); | |
3bb122c7 | 157 | hMult->Reset(); |
158 | } | |
159 | ||
160 | } | |
161 | ||
162 | ||
163 | for(UShort_t det=1;det<=3;det++) { | |
3bb122c7 | 164 | Int_t nRings = (det==1 ? 1 : 2); |
165 | for (UShort_t ir = 0; ir < nRings; ir++) { | |
bb8a464f | 166 | |
3bb122c7 | 167 | Char_t ring = (ir == 0 ? 'I' : 'O'); |
9f55be54 | 168 | TH2F* hMult = (TH2F*)fOutputList->FindObject(Form("FMD%d%c_vtxbin%d",det,ring,vtxbin)); |
169 | ||
3bb122c7 | 170 | UShort_t nsec = (ir == 0 ? 20 : 40); |
171 | UShort_t nstr = (ir == 0 ? 512 : 256); | |
da0805e2 | 172 | /* |
173 | TH1F* hEnergyDist = pars->GetEmptyEnergyDistribution(det,ring); | |
174 | TF1* fitFunc = hEnergyDist->GetFunction("FMDfitFunc"); | |
175 | TH1F* hSignalDist = pars->GetRingEnergyDistribution(det, ring); | |
176 | TF1* fitFuncSignal = hSignalDist->GetFunction("FMDfitFunc"); | |
177 | */ | |
3bb122c7 | 178 | for(UShort_t sec =0; sec < nsec; sec++) { |
179 | for(UShort_t strip = 0; strip < nstr; strip++) { | |
7c3e5162 | 180 | Float_t mult = fESD->Multiplicity(det,ring,sec,strip); |
9f55be54 | 181 | |
bb8a464f | 182 | if(mult == 0 || mult == AliESDFMD::kInvalidMult) continue; |
9f55be54 | 183 | |
f58a4769 | 184 | Float_t phi = pars->GetPhiFromSector(det,ring,sec); |
185 | Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex[2]); | |
41bad769 | 186 | |
da0805e2 | 187 | Float_t mult_cut = 0.3;//pars->GetMPV(det,ring,eta);// - pars->GetSigma(det,ring,eta);//0.25;//0.15;//m-2*s;//0.15;//0.2;//m-3*s;// 0.2;//0.01;//m-2*s;//0.2; |
36f7d459 | 188 | // if(ring == 'I') |
189 | // mult_cut = 0.10; | |
b85ea106 | 190 | //Float_t mult_cut = pars->GetMPV(det,ring,eta) - 5*pars->GetSigma(det,ring,eta); |
bb8a464f | 191 | Float_t nParticles = 0; |
6289b3e8 | 192 | if(fESD->GetUniqueID() == kTRUE) { |
193 | //proton + proton | |
cc066cb9 | 194 | |
195 | if(mult > mult_cut) { | |
da0805e2 | 196 | nParticles = 1.; |
cc066cb9 | 197 | } |
6289b3e8 | 198 | } |
199 | else { | |
200 | ||
201 | //Pb+Pb | |
202 | Float_t mpv = pars->GetMPV(det,ring,eta); | |
203 | Float_t sigma = pars->GetSigma(det,ring,eta); | |
204 | Float_t alpha = pars->Get2MIPWeight(det,ring,eta); | |
205 | Float_t beta = pars->Get3MIPWeight(det,ring,eta); | |
206 | ||
207 | Float_t sumCor = TMath::Landau(mult,mpv,sigma,kTRUE)+ | |
208 | alpha*TMath::Landau(mult,2*mpv+2*sigma*TMath::Log(2),2*sigma,kTRUE)+ | |
209 | beta*TMath::Landau(mult,3*mpv+3*sigma*TMath::Log(3),3*sigma,kTRUE); | |
210 | Float_t weight = TMath::Landau(mult,mpv,sigma,kTRUE)+ | |
211 | 2*alpha*TMath::Landau(mult,2*mpv+2*sigma*TMath::Log(2),2*sigma,kTRUE)+ | |
212 | 3*beta*TMath::Landau(mult,3*mpv+3*sigma*TMath::Log(3),3*sigma,kTRUE); | |
213 | ||
214 | ||
cc066cb9 | 215 | if(mult > mult_cut) { |
6289b3e8 | 216 | if(sumCor) nParticles = weight / sumCor; |
217 | else nParticles = 1; | |
cc066cb9 | 218 | |
6289b3e8 | 219 | } |
220 | //std::cout<<sumCor<<" "<<weight<<" "<<" "<<mult<<" "<<nParticles<<std::endl; | |
221 | ||
222 | } | |
223 | ||
224 | ||
cc066cb9 | 225 | |
226 | ||
227 | ||
1282ce49 | 228 | Float_t correction = GetAcceptanceCorrection(ring,strip); |
f58a4769 | 229 | |
230 | //std::cout<<"before "<<correction<<std::endl; | |
78f6f750 | 231 | if(fESD->GetUniqueID() == kTRUE) { |
4fb49e43 | 232 | TH1F* hDoubleHitCorrection = pars->GetDoubleHitCorrection(det,ring); |
41bad769 | 233 | |
4fb49e43 | 234 | if(hDoubleHitCorrection->GetBinContent(hDoubleHitCorrection->FindBin(eta)) != 0) |
b85ea106 | 235 | correction = correction*hDoubleHitCorrection->GetBinContent(hDoubleHitCorrection->FindBin(eta)); |
41bad769 | 236 | |
78f6f750 | 237 | } |
f58a4769 | 238 | |
f7356393 | 239 | //Dead channel correction: |
240 | TH1F* hFMDDeadCorrection = pars->GetFMDDeadCorrection(vtxbin); | |
241 | if(hFMDDeadCorrection) | |
242 | if(hFMDDeadCorrection->GetBinContent(hFMDDeadCorrection->FindBin(eta)) != 0) | |
243 | correction = correction*hFMDDeadCorrection->GetBinContent(hFMDDeadCorrection->FindBin(eta)); | |
244 | ||
245 | ||
da0805e2 | 246 | //std::cout<<det<<" "<<ring<<" "<<sec<<" "<<strip<<" "<<vertex[2]<<" "<<eta<<std::endl; |
247 | ||
248 | //Float_t signal = hSignalDist->Integral(hSignalDist->FindBin(0.5),hSignalDist->FindBin(2)) ;//pars->GetConstant(det,ring,eta); | |
249 | //Float_t bkg = hEnergyDist->GetBinContent(hEnergyDist->FindBin(mult)); | |
250 | //Float_t signal = hSignalDist->GetBinContent(hSignalDist->FindBin(mult)); | |
251 | /* | |
252 | if(fitFunc && fitFuncSignal && pars->IsRealData()) { | |
253 | Float_t bkg = fitFunc->Eval(mult); | |
254 | Float_t signal = fitFuncSignal->Eval(mult); | |
255 | ||
256 | if(signal > 0) { | |
257 | correction = correction/(1-(bkg/signal)); | |
258 | //test | |
259 | // if(TMath::Abs(eta)<3) correction = 1.15*correction; | |
260 | // if(det == 2 && ring == 'I') correction = correction/(1-bkg/signal); | |
261 | //if(det == 2 && ring == 'O') correction = correction/(1-bkg/signal); | |
262 | //if(det == 3 && ring == 'I') correction = correction/(1-bkg/signal); | |
263 | //if(det == 3 && ring == 'O') correction = correction/(1-bkg/signal); | |
264 | ||
265 | } | |
266 | }*/ | |
6289b3e8 | 267 | if(correction) nParticles = nParticles / correction; |
cc066cb9 | 268 | if(nParticles > 0) |
269 | hMult->Fill(eta,phi,nParticles); | |
bb8a464f | 270 | |
3bb122c7 | 271 | |
272 | } | |
273 | } | |
bb8a464f | 274 | |
3bb122c7 | 275 | } |
bb8a464f | 276 | |
3bb122c7 | 277 | |
278 | ||
279 | } | |
78f6f750 | 280 | |
da0805e2 | 281 | |
78f6f750 | 282 | |
7c3e5162 | 283 | if(fStandalone) { |
284 | PostData(0, fOutputList); | |
285 | } | |
3bb122c7 | 286 | |
287 | } | |
288 | //_____________________________________________________________________ | |
1282ce49 | 289 | Float_t AliFMDAnalysisTaskDensity::GetAcceptanceCorrection(Char_t ring, UShort_t strip) |
290 | { | |
78f6f750 | 291 | AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance(); |
1282ce49 | 292 | |
78f6f750 | 293 | //AliFMDRing fmdring(ring); |
294 | //fmdring.Init(); | |
295 | Float_t rad = pars->GetMaxR(ring)-pars->GetMinR(ring); | |
296 | Float_t nstrips = (ring == 'I' ? 512 : 256); | |
297 | Float_t segment = rad / nstrips; | |
298 | Float_t radius = pars->GetMinR(ring) + segment*strip; | |
299 | ||
300 | Float_t basearea1 = 0.5*pars->GetBaseStripLength(ring,strip)*TMath::Power(radius,2); | |
301 | Float_t basearea2 = 0.5*pars->GetBaseStripLength(ring,strip)*TMath::Power((radius-segment),2); | |
1282ce49 | 302 | Float_t basearea = basearea1 - basearea2; |
303 | ||
78f6f750 | 304 | Float_t area1 = 0.5*pars->GetStripLength(ring,strip)*TMath::Power(radius,2); |
305 | Float_t area2 = 0.5*pars->GetStripLength(ring,strip)*TMath::Power((radius-segment),2); | |
1282ce49 | 306 | Float_t area = area1 - area2; |
307 | ||
308 | Float_t correction = area/basearea; | |
309 | ||
310 | return correction; | |
311 | } | |
f58a4769 | 312 | //_____________________________________________________________________ |
3bb122c7 | 313 | // |
314 | //EOF | |
315 | // |