]>
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 "TH1F.h" | |
10 | #include "TH2F.h" | |
11 | #include "AliFMDAnalysisTaskDndeta.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 "AliLog.h" | |
20 | #include "AliESDVertex.h" | |
21 | #include "TMath.h" | |
22 | #include "AliFMDAnaParameters.h" | |
23 | //#include "AliFMDGeometry.h" | |
24 | #include "AliGenEventHeader.h" | |
25 | #include "AliGenPythiaEventHeader.h" | |
26 | #include "AliHeader.h" | |
27 | //#include "TDatabasePDG.h" | |
28 | //#include "TParticlePDG.h" | |
29 | #include "AliFMDStripIndex.h" | |
30 | #include "AliESDInputHandler.h" | |
31 | ClassImp(AliFMDAnalysisTaskDndeta) | |
32 | ||
33 | ||
34 | AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta() | |
35 | : fDebug(0), | |
36 | fOutputList(0), | |
37 | fInputList(0), | |
38 | fVertexString(0x0), | |
39 | fNevents(), | |
40 | fNMCevents(), | |
41 | fStandalone(kTRUE), | |
42 | fLastTrackByStrip(0) | |
43 | { | |
44 | // Default constructor | |
45 | DefineInput (0, TList::Class()); | |
46 | DefineOutput(0, TList::Class()); | |
47 | } | |
48 | //_____________________________________________________________________ | |
49 | AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta(const char* name, Bool_t SE): | |
50 | AliAnalysisTask(name, "Density"), | |
51 | fDebug(0), | |
52 | fOutputList(0), | |
53 | fInputList(0), | |
54 | fVertexString(0x0), | |
55 | fNevents(), | |
56 | fNMCevents(), | |
57 | fStandalone(kTRUE), | |
58 | fLastTrackByStrip(0) | |
59 | { | |
60 | fStandalone = SE; | |
61 | if(fStandalone) { | |
62 | DefineInput (0, TList::Class()); | |
63 | DefineInput(1, TObjString::Class()); | |
64 | DefineOutput(0, TList::Class()); | |
65 | ||
66 | } | |
67 | } | |
68 | //_____________________________________________________________________ | |
69 | void AliFMDAnalysisTaskDndeta::CreateOutputObjects() | |
70 | { | |
71 | AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance(); | |
72 | ||
73 | if(!fOutputList) | |
74 | fOutputList = new TList(); | |
75 | fOutputList->SetName("BackgroundCorrected"); | |
76 | ||
77 | ||
78 | TH2F* hMult = 0; | |
79 | TH1F* hHits = 0; | |
80 | TH2F* hMultTrVtx = 0; | |
81 | TH1F* hPrimVertexBin = 0; | |
82 | ||
83 | ||
84 | TH2F* hBgTmp = pars->GetBackgroundCorrection(1, 'I', 0); | |
85 | TH1F* hPrimary = new TH1F("hMultvsEta","hMultvsEta", | |
86 | hBgTmp->GetNbinsX(), | |
87 | hBgTmp->GetXaxis()->GetXmin(), | |
88 | hBgTmp->GetXaxis()->GetXmax()); | |
89 | hPrimary->Sumw2(); | |
90 | fOutputList->Add(hPrimary); | |
91 | Int_t nVtxbins = pars->GetNvtxBins(); | |
92 | TH2F* hBg = 0; | |
93 | for(Int_t i = 0; i< nVtxbins; i++) { | |
94 | ||
95 | for(Int_t det =1; det<=3;det++) | |
96 | { | |
97 | Int_t nRings = (det==1 ? 1 : 2); | |
98 | for(Int_t ring = 0;ring<nRings;ring++) | |
99 | { | |
100 | Char_t ringChar = (ring == 0 ? 'I' : 'O'); | |
101 | Int_t nSec = (ring == 0 ? 20 : 40); | |
102 | ||
103 | ||
104 | ||
105 | hBg = pars->GetBackgroundCorrection(det, ringChar, i); | |
106 | hMult = new TH2F(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),Form("dNdeta_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()); | |
111 | hMultTrVtx = new TH2F(Form("dNdetaTrVtx_FMD%d%c_vtxbin%d",det,ringChar,i),Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i), | |
112 | hBg->GetNbinsX(), | |
113 | hBg->GetXaxis()->GetXmin(), | |
114 | hBg->GetXaxis()->GetXmax(), | |
115 | nSec, 0, 2*TMath::Pi()); | |
116 | hHits = new TH1F(Form("hMCHits_FMD%d%c_vtxbin%d",det,ringChar,i),Form("hMCHits_FMD%d%c_vtxbin%d",det,ringChar,i), | |
117 | hBg->GetNbinsX(), | |
118 | hBg->GetXaxis()->GetXmin(), | |
119 | hBg->GetXaxis()->GetXmax()); | |
120 | hHits->Sumw2(); | |
121 | fOutputList->Add(hHits); | |
122 | ||
123 | hMult->Sumw2(); | |
124 | fOutputList->Add(hMult); | |
125 | ||
126 | hMultTrVtx->Sumw2(); | |
127 | fOutputList->Add(hMultTrVtx); | |
128 | ||
129 | } | |
130 | } | |
131 | } | |
132 | ||
133 | for(Int_t i = 0; i< nVtxbins; i++) { | |
134 | ||
135 | hPrimVertexBin = new TH1F(Form("primmult_vtxbin%d",i), | |
136 | Form("primmult_vtxbin%d",i), | |
137 | hBgTmp->GetNbinsX(), | |
138 | hBgTmp->GetXaxis()->GetXmin(), | |
139 | hBgTmp->GetXaxis()->GetXmax()); | |
140 | hPrimVertexBin->Sumw2(); | |
141 | fOutputList->Add(hPrimVertexBin); | |
142 | //SPD part | |
143 | TH2F* hSPDMult = new TH2F(Form("dNdeta_SPD_vtxbin%d",i),Form("dNdeta_SPD_vtxbin%d",i), | |
144 | hBgTmp->GetNbinsX(), | |
145 | hBgTmp->GetXaxis()->GetXmin(), | |
146 | hBgTmp->GetXaxis()->GetXmax(), | |
147 | 20, 0, 2*TMath::Pi()); | |
148 | hSPDMult->Sumw2(); | |
149 | fOutputList->Add(hSPDMult); | |
150 | TH2F* hSPDMultTrVtx = new TH2F(Form("dNdetaTrVtx_SPD_vtxbin%d",i),Form("dNdetaTrVtx_SPD_vtxbin%d",i), | |
151 | hBgTmp->GetNbinsX(), | |
152 | hBgTmp->GetXaxis()->GetXmin(), | |
153 | hBgTmp->GetXaxis()->GetXmax(), | |
154 | 20, 0, 2*TMath::Pi()); | |
155 | hSPDMultTrVtx->Sumw2(); | |
156 | fOutputList->Add(hSPDMultTrVtx); | |
157 | } | |
158 | ||
159 | fNevents.SetBins(nVtxbins,0,nVtxbins); | |
160 | fNevents.SetName("nEvents"); | |
161 | fNMCevents.SetBins(nVtxbins,0,nVtxbins); | |
162 | fNMCevents.SetName("nMCEvents"); | |
163 | fOutputList->Add(&fNevents); | |
164 | fOutputList->Add(&fNMCevents); | |
165 | ||
166 | } | |
167 | //_____________________________________________________________________ | |
168 | void AliFMDAnalysisTaskDndeta::ConnectInputData(Option_t */*option*/) | |
169 | { | |
170 | if(fStandalone) { | |
171 | fInputList = (TList*)GetInputData(0); | |
172 | } | |
173 | } | |
174 | //_____________________________________________________________________ | |
175 | void AliFMDAnalysisTaskDndeta::Exec(Option_t */*option*/) | |
176 | { | |
177 | AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance(); | |
178 | fVertexString = (TObjString*)fInputList->At(0); | |
179 | ||
180 | Int_t vtxbin = fVertexString->GetString().Atoi(); | |
181 | ||
182 | fNevents.Fill(vtxbin); | |
183 | ||
184 | for(UShort_t det=1;det<=3;det++) { | |
185 | Int_t nRings = (det==1 ? 1 : 2); | |
186 | for (UShort_t ir = 0; ir < nRings; ir++) { | |
187 | Char_t ringChar = (ir == 0 ? 'I' : 'O'); | |
188 | ||
189 | TH2F* hMultTotal = (TH2F*)fOutputList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,vtxbin)); | |
190 | TH2F* hMultTotalTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_FMD%d%c_vtxbin%d",det,ringChar,vtxbin)); | |
191 | ||
192 | TH2F* hMultInput = (TH2F*)fInputList->FindObject(Form("mult_FMD%d%c_vtxbin%d",det,ringChar,vtxbin)); | |
193 | TH2F* hMultInputTrVtx = (TH2F*)fInputList->FindObject(Form("multTrVtx_FMD%d%c_vtxbin%d",det,ringChar,vtxbin)); | |
194 | ||
195 | hMultTotal->Add(hMultInput); | |
196 | hMultTotalTrVtx->Add(hMultInputTrVtx); | |
197 | ||
198 | } | |
199 | } | |
200 | ||
201 | //SPD code | |
202 | TH2F* hMultSPDTotal = (TH2F*)fOutputList->FindObject(Form("dNdeta_SPD_vtxbin%d",vtxbin)); | |
203 | TH2F* hMultSPDTotalTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_SPD_vtxbin%d",vtxbin)); | |
204 | TH2F* hMultSPDInput = (TH2F*)fInputList->FindObject(Form("mult_SPD_vtxbin%d",vtxbin)); | |
205 | TH2F* hMultSPDInputTrVtx = (TH2F*)fInputList->FindObject(Form("multTrVtx_SPD_vtxbin%d",vtxbin)); | |
206 | hMultSPDTotal->Add(hMultSPDInput); | |
207 | hMultSPDTotalTrVtx->Add(hMultSPDInputTrVtx); | |
208 | ||
209 | ||
210 | if(pars->GetProcessPrimary()) | |
211 | ProcessPrimary(); | |
212 | ||
213 | if(fStandalone) { | |
214 | PostData(0, fOutputList); | |
215 | } | |
216 | ||
217 | } | |
218 | //_____________________________________________________________________ | |
219 | void AliFMDAnalysisTaskDndeta::Terminate(Option_t */*option*/) { | |
220 | ||
221 | AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance(); | |
222 | ||
223 | Int_t nVtxbins = pars->GetNvtxBins(); | |
224 | ||
225 | for(UShort_t det=1;det<=3;det++) { | |
226 | Int_t nRings = (det==1 ? 1 : 2); | |
227 | for (UShort_t ir = 0; ir < nRings; ir++) { | |
228 | Char_t ringChar = (ir == 0 ? 'I' : 'O'); | |
229 | for(Int_t i =0; i<nVtxbins; i++) { | |
230 | ||
231 | TH2F* hMultTotal = (TH2F*)fOutputList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i)); | |
232 | //TH2F* hMultTrVtx = (TH2F*)hMultTotal->Clone(Form("dNdeta_FMD%d%c_TrVtx_vtxbin%d",det,ringChar,i)); | |
233 | TH2F* hMultTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_FMD%d%c_vtxbin%d",det,ringChar,i)); | |
234 | ||
235 | TH1D* hMultProj = hMultTotal->ProjectionX(Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ringChar,i),1,hMultTotal->GetNbinsY()); | |
236 | TH1D* hMultProjTrVtx = hMultTrVtx->ProjectionX(Form("dNdeta_FMD%d%c_TrVtx_vtxbin%d_proj",det,ringChar,i),1,hMultTotal->GetNbinsY()); | |
237 | //fOutputList->Add(hMultTrVtx); | |
238 | fOutputList->Add(hMultProj); | |
239 | fOutputList->Add(hMultProjTrVtx); | |
240 | } | |
241 | } | |
242 | } | |
243 | ||
244 | for(Int_t i =0; i<nVtxbins; i++) { | |
245 | ||
246 | TH2F* hSPDMult = (TH2F*)fOutputList->FindObject(Form("dNdeta_SPD_vtxbin%d",i)); | |
247 | TH2F* hSPDMultTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_SPD_vtxbin%d",i)); | |
248 | ||
249 | TH1D* hMultProj = hSPDMult->ProjectionX(Form("dNdeta_SPD_vtxbin%d_proj",i),1,hSPDMult->GetNbinsY()); | |
250 | TH1D* hMultProjTrVtx = hSPDMultTrVtx->ProjectionX(Form("dNdetaTrVtx_SPD_vtxbin%d_proj",i),1,hSPDMultTrVtx->GetNbinsY()); | |
251 | ||
252 | fOutputList->Add(hMultProj); | |
253 | fOutputList->Add(hMultProjTrVtx); | |
254 | ||
255 | } | |
256 | ||
257 | std::cout<<"FMD analysis accepted "<<fNevents.GetEntries()<<" events"<<std::endl; | |
258 | } | |
259 | //_____________________________________________________________________ | |
260 | void AliFMDAnalysisTaskDndeta::ProcessPrimary() { | |
261 | ||
262 | AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler()); | |
263 | AliMCEvent* mcEvent = eventHandler->MCEvent(); | |
264 | if(!mcEvent) | |
265 | return; | |
266 | ||
267 | fLastTrackByStrip.Reset(-1); | |
268 | ||
269 | AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance(); | |
270 | ||
271 | AliMCParticle* particle = 0; | |
272 | AliStack* stack = mcEvent->Stack(); | |
273 | ||
274 | TH1F* hPrimary = (TH1F*)fOutputList->FindObject("hMultvsEta"); | |
275 | AliHeader* header = mcEvent->Header(); | |
276 | AliGenEventHeader* genHeader = header->GenEventHeader(); | |
277 | ||
278 | TArrayF vertex; | |
279 | genHeader->PrimaryVertex(vertex); | |
280 | if(TMath::Abs(vertex.At(2)) > pars->GetVtxCutZ()) | |
281 | return; | |
282 | Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins(); | |
283 | Double_t vertexBinDouble = (vertex.At(2) + pars->GetVtxCutZ()) / delta; | |
284 | Int_t vertexBin = (Int_t)vertexBinDouble; | |
285 | ||
286 | Bool_t firstTrack = kTRUE; | |
287 | ||
288 | // we loop over the primaries only unless we need the hits (diagnostics running slowly) | |
289 | Int_t nTracks = stack->GetNprimary(); | |
290 | if(pars->GetProcessHits()) | |
291 | nTracks = stack->GetNtrack(); | |
292 | ||
293 | for(Int_t i = 0 ;i<nTracks;i++) { | |
294 | particle = (AliMCParticle*) mcEvent->GetTrack(i); | |
295 | if(!particle) | |
296 | continue; | |
297 | ||
298 | if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) { | |
299 | hPrimary->Fill(particle->Eta()); | |
300 | ||
301 | ||
302 | TH1F* hPrimVtxBin = (TH1F*)fOutputList->FindObject(Form("primmult_vtxbin%d",vertexBin)); | |
303 | hPrimVtxBin->Fill(particle->Eta()); | |
304 | if(firstTrack) { | |
305 | fNMCevents.Fill(vertexBin); | |
306 | firstTrack = kFALSE; | |
307 | } | |
308 | ||
309 | } | |
310 | if(pars->GetProcessHits()) { | |
311 | ||
312 | for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) { | |
313 | ||
314 | AliTrackReference* ref = particle->GetTrackReference(j); | |
315 | UShort_t det,sec,strip; | |
316 | Char_t ring; | |
317 | if(ref->DetectorId() != AliTrackReference::kFMD) | |
318 | continue; | |
319 | AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip); | |
320 | Float_t thisStripTrack = fLastTrackByStrip(det,ring,sec,strip); | |
321 | if(particle->Charge() != 0 && i != thisStripTrack ) { | |
322 | //Double_t x,y,z; | |
323 | ||
324 | Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));//-1*TMath::Log(TMath::Tan(0.5*theta)); | |
325 | TH1F* hHits = (TH1F*)fOutputList->FindObject(Form("hMCHits_FMD%d%c_vtxbin%d",det,ring,vertexBin)); | |
326 | ||
327 | ||
328 | hHits->Fill(eta); | |
329 | ||
330 | Float_t nstrips = (ring =='O' ? 256 : 512); | |
331 | ||
332 | fLastTrackByStrip(det,ring,sec,strip) = (Float_t)i; | |
333 | ||
334 | if(strip >0) | |
335 | fLastTrackByStrip(det,ring,sec,strip-1) = (Float_t)i; | |
336 | if(strip < (nstrips - 1)) | |
337 | fLastTrackByStrip(det,ring,sec,strip+1) = (Float_t)i; | |
338 | ||
339 | } | |
340 | ||
341 | ||
342 | } | |
343 | ||
344 | ||
345 | } | |
346 | ||
347 | ||
348 | ||
349 | } | |
350 | } | |
351 | //_____________________________________________________________________ | |
352 | // | |
353 | // | |
354 | // EOF |