Next round of upgrades and cleanups of code. The code is now more streamlined and...
[u/mrichter/AliRoot.git] / FMD / analysis / AliFMDAnalysisTaskGenerateCorrection.cxx
CommitLineData
d31dce33 1#include "AliFMDAnalysisTaskGenerateCorrection.h"
2#include "AliESDEvent.h"
3#include "iostream"
4#include "AliESDFMD.h"
5#include "TH2F.h"
6#include "AliTrackReference.h"
7#include "AliStack.h"
8#include "AliFMDAnaParameters.h"
9#include "AliFMDStripIndex.h"
10#include "AliStack.h"
11#include "AliMCParticle.h"
12#include "AliMCEvent.h"
13//#include "AliFMDGeometry.h"
14#include "TArray.h"
15#include "AliGenEventHeader.h"
16#include "AliHeader.h"
17#include "AliFMDAnaCalibBackgroundCorrection.h"
18#include "AliFMDAnaCalibEventSelectionEfficiency.h"
19//#include "AliCDBManager.h"
20//#include "AliCDBId.h"
21//#include "AliCDBMetaData.h"
22#include "TSystem.h"
23#include "TROOT.h"
24#include "TAxis.h"
25ClassImp(AliFMDAnalysisTaskGenerateCorrection)
26
27//_____________________________________________________________________
28AliFMDAnalysisTaskGenerateCorrection::AliFMDAnalysisTaskGenerateCorrection():
29AliAnalysisTaskSE(),
30 fListOfHits(),
31 fListOfPrimaries(),
32 fListOfCorrection(),
33 fVertexBins(),
34 fLastTrackByStrip(0),
35 fHitsByStrip(0),
36 fZvtxCut(10),
37 fNvtxBins(10),
38 fNbinsEta(200),
9f55be54 39 fBackground(0),
40 fEventSelectionEff(0)
d31dce33 41{
42 // Default constructor
43}
44//_____________________________________________________________________
45AliFMDAnalysisTaskGenerateCorrection::AliFMDAnalysisTaskGenerateCorrection(const char* name):
46 AliAnalysisTaskSE(name),
47 fListOfHits(),
48 fListOfPrimaries(),
49 fListOfCorrection(),
50 fVertexBins(),
51 fLastTrackByStrip(0),
52 fHitsByStrip(0),
53 fZvtxCut(10),
54 fNvtxBins(10),
55 fNbinsEta(200),
9f55be54 56 fBackground(0),
57 fEventSelectionEff(0)
d31dce33 58{
59
60 DefineOutput(1, TList::Class());
61 DefineOutput(2, TList::Class());
62 DefineOutput(3, TH1F::Class());
63 DefineOutput(4, TList::Class());
64}
65//_____________________________________________________________________
66void AliFMDAnalysisTaskGenerateCorrection::UserCreateOutputObjects()
67{
68// Create the output containers
69//
70
71 std::cout<<"Creating output objects"<<std::endl;
72 for(Int_t iring = 0; iring<2;iring++) {
73 Char_t ringChar = (iring == 0 ? 'I' : 'O');
74 Int_t nSec = (iring == 1 ? 40 : 20);
75 for(Int_t v=0; v<fNvtxBins;v++) {
76
77 TH2F* hPrimary = new TH2F(Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
78 Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
79 fNbinsEta, -6,6, nSec, 0,2*TMath::Pi());
80 hPrimary->Sumw2();
81 fListOfPrimaries.Add(hPrimary);
82 }
83 }
84
85
86
87 for(Int_t det =1; det<=3;det++) {
88 Int_t nRings = (det==1 ? 1 : 2);
89 for(Int_t ring = 0;ring<nRings;ring++) {
90 Int_t nSec = (ring == 1 ? 40 : 20);
91 Char_t ringChar = (ring == 0 ? 'I' : 'O');
92 TH1F* doubleHits = new TH1F(Form("DoubleHits_FMD%d%c",det,ringChar),
93 Form("DoubleHits_FMD%d%c",det,ringChar),fNbinsEta, -6,6);
94 TH1F* allHits = new TH1F(Form("allHits_FMD%d%c",det,ringChar),
95 Form("allHits_FMD%d%c",det,ringChar), fNbinsEta, -6,6);
96
97 doubleHits->Sumw2();
98 allHits->Sumw2();
99 fListOfHits.Add(allHits);
100 fListOfHits.Add(doubleHits);
101
102 for(Int_t v=0; v<fNvtxBins;v++) {
103 TH2F* hHits = new TH2F(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
104 Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
105 fNbinsEta, -6,6, nSec, 0, 2*TMath::Pi());
106 hHits->Sumw2();
107 fListOfHits.Add(hHits);
108
109 }
110 }
111 }
112
113 TH1F* hEventsSelected = new TH1F("EventsSelected","EventsSelected",fNvtxBins,0,fNvtxBins);
114 TH1F* hEventsAll = new TH1F("EventsAll","EventsAll",fNvtxBins,0,fNvtxBins);
115 // TH1F* hTriggered = new TH1F("Triggered","Triggered",fNvtxBins,0,fNvtxBins);
116 // TH1F* hTriggeredAll = new TH1F("TriggeredAll","TriggeredAll",fNvtxBins,0,fNvtxBins);
117
118 fListOfHits.Add(hEventsSelected);
119 fListOfPrimaries.Add(hEventsAll);
120
121 // fListOfHits.Add(hTriggered);
122 // fListOfPrimaries.Add(hTriggeredAll);
123
124 fVertexBins.SetName("VertexBins");
125 fVertexBins.GetXaxis()->Set(fNvtxBins,-1*fZvtxCut,fZvtxCut);
126
127}
128//_____________________________________________________________________
129void AliFMDAnalysisTaskGenerateCorrection::Init()
130{
131 fLastTrackByStrip.Reset(-1);
132
133
134}
135//_____________________________________________________________________
136void AliFMDAnalysisTaskGenerateCorrection::UserExec(Option_t */*option*/)
137{
138
139 fLastTrackByStrip.Reset(-1);
140 fHitsByStrip.Reset(0);
141 AliMCEvent* mcevent = MCEvent();
142
143 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
144
145
146 AliESDEvent* esdevent = (AliESDEvent*)InputEvent();
147 Double_t esdvertex[3];
148 pars->GetVertex(esdevent,esdvertex);
149
150 AliMCParticle* particle = 0;
151 AliStack* stack = mcevent->Stack();
152
153 UShort_t det,sec,strip;
154 Char_t ring;
155
156 Int_t nTracks = mcevent->GetNumberOfTracks();
157 AliHeader* header = mcevent->Header();
158 AliGenEventHeader* genHeader = header->GenEventHeader();
159
160 TArrayF vertex;
161 genHeader->PrimaryVertex(vertex);
162
163 if(TMath::Abs(vertex.At(2)) > fZvtxCut)
164 return;
165
166 Double_t delta = 2*fZvtxCut/fNvtxBins;
167 Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
168 Int_t vertexBin = (Int_t)vertexBinDouble;
169
170 // Vertex determination correction
171 TH1F* hEventsSelected = (TH1F*)fListOfHits.FindObject("EventsSelected");
172 TH1F* hEventsAll = (TH1F*)fListOfPrimaries.FindObject("EventsAll");
173 // TH1F* hTriggered = (TH1F*)fListOfHits.FindObject("Triggered");
174 // TH1F* hTriggeredAll = (TH1F*)fListOfPrimaries.FindObject("TriggeredAll");
175
176 Bool_t vtxFound = kTRUE;
177 if(esdvertex[0] == 0 && esdvertex[1] == 0 && esdvertex[2] == 0)
178 vtxFound = kFALSE;
179
180 Bool_t isTriggered = pars->IsEventTriggered(esdevent);
181
182 if(vtxFound && isTriggered) {
183 //hTriggered->Fill(vertexBin);
184 hEventsSelected->Fill(vertexBin);
185 }
186
187 //if(isTriggered)
188
189 hEventsAll->Fill(vertexBin);
190 //if(vtxFound)
191 // hTriggeredAll->Fill(vertexBin);
192
193 for(Int_t i = 0 ;i<nTracks;i++) {
194 particle = mcevent->GetTrack(i);
195
196 if(!particle)
197 continue;
198
199 if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
200
201
202 TH2F* hPrimaryInner = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",'I',vertexBin));
203 TH2F* hPrimaryOuter = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",'O',vertexBin));
204 hPrimaryInner->Fill(particle->Eta(),particle->Phi());
205 hPrimaryOuter->Fill(particle->Eta(),particle->Phi());
206 }
207
208 for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) {
209
210 AliTrackReference* ref = particle->GetTrackReference(j);
211
212 if(ref->DetectorId() != AliTrackReference::kFMD)
213 continue;
214 AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip);
215 Float_t thisStripTrack = fLastTrackByStrip(det,ring,sec,strip);
216 if(particle->Charge() != 0 && i != thisStripTrack ) {
217
218 Float_t phi = pars->GetPhiFromSector(det,ring,sec);
219 Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));
220
221 TH2F* hHits = (TH2F*)fListOfHits.FindObject(Form("hHits_FMD%d%c_vtx%d", det,ring,vertexBin));
222 hHits->Fill(eta,phi);
223 Float_t nstrips = (ring =='O' ? 256 : 512);
224 fHitsByStrip(det,ring,sec,strip) +=1;
225 TH1F* allHits = (TH1F*)fListOfHits.FindObject(Form("allHits_FMD%d%c",det,ring));
226 TH1F* doubleHits = (TH1F*)fListOfHits.FindObject(Form("DoubleHits_FMD%d%c",det,ring));
227
228 if(fHitsByStrip(det,ring,sec,strip) == 1)
229 allHits->Fill(eta);
230
231 doubleHits->Fill(eta);
232
233 fLastTrackByStrip(det,ring,sec,strip) = (Float_t)i;
234 if(strip >0)
235 fLastTrackByStrip(det,ring,sec,strip-1) = (Float_t)i;
236 if(strip < (nstrips - 1))
237 fLastTrackByStrip(det,ring,sec,strip+1) = (Float_t)i;
238 }
239 }
240
241 }
242
243
244 PostData(1, &fListOfHits);
245 PostData(2, &fListOfPrimaries);
246 PostData(3, &fVertexBins);
247}
248//_____________________________________________________________________
249void AliFMDAnalysisTaskGenerateCorrection::Terminate(Option_t */*option*/)
250{
251 /* TH1F* allHits = (TH1F*)fListOfHits.FindObject("allHits");
252 TH1F* doubleHits = (TH1F*)fListOfHits.FindObject("DoubleHits");
253
254 doubleHits->Divide(allHits);
255 GenerateCorrection();
256 PostData(1, &fListOfHits);
257 PostData(4, &fListOfCorrection);*/
258
259}
260//_____________________________________________________________________
261void AliFMDAnalysisTaskGenerateCorrection::GenerateCorrection() {
262
263 fBackground = new AliFMDAnaCalibBackgroundCorrection();
264 fEventSelectionEff = new AliFMDAnaCalibEventSelectionEfficiency();
265
266 //TH1F* hTriggered = (TH1F*)fListOfHits.FindObject("Triggered");
267 //TH1F* hTriggeredAll = (TH1F*)fListOfPrimaries.FindObject("TriggeredAll");
268
269 TH1F* hEventsSelected = (TH1F*)fListOfHits.FindObject("EventsSelected");
270 TH1F* hEventsAll = (TH1F*)fListOfPrimaries.FindObject("EventsAll");
271
272 hEventsSelected->Divide(hEventsAll);
273
274 fEventSelectionEff->SetCorrection(hEventsSelected);
275
276 for(Int_t det= 1; det <=3; det++) {
277 Int_t nRings = (det==1 ? 1 : 2);
278
279 for(Int_t iring = 0; iring<nRings; iring++) {
280 Char_t ring = (iring == 0 ? 'I' : 'O');
281 TH1F* allHits = (TH1F*)fListOfHits.FindObject(Form("allHits_FMD%d%c",det,ring));
282 TH1F* doubleHits = (TH1F*)fListOfHits.FindObject(Form("DoubleHits_FMD%d%c",det,ring));
283 fBackground->SetDoubleHitCorrection(det,ring,doubleHits);
284 doubleHits->Divide(allHits);
285 for(Int_t vertexBin=0;vertexBin<fNvtxBins ;vertexBin++) {
286 TH2F* hHits = (TH2F*)fListOfHits.FindObject(Form("hHits_FMD%d%c_vtx%d", det,ring,vertexBin));
287 TH2F* hPrimary = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",ring,vertexBin));
288 TH2F* hCorrection = (TH2F*)hHits->Clone(Form("FMD%d%c_vtxbin_%d_correction",det,ring,vertexBin));
289 hCorrection->Divide(hPrimary);
290 hCorrection->SetTitle(hCorrection->GetName());
291 fListOfCorrection.Add(hCorrection);
292 fBackground->SetBgCorrection(det,ring,vertexBin,hCorrection);
293
294
295 }
296
297 }
298 }
299 TAxis refAxis(fNvtxBins,-1*fZvtxCut,fZvtxCut);
300 fBackground->SetRefAxis(&refAxis);
301
302}
303//_____________________________________________________________________
304void AliFMDAnalysisTaskGenerateCorrection::ReadFromFile(const Char_t* filename, Bool_t storeInOCDB, Int_t /*runNo*/) {
305
306 TFile infile(filename);
307 TH1F* hVertex = (TH1F*)infile.Get("VertexBins");
308 fZvtxCut = hVertex->GetXaxis()->GetXmax();
309 fNvtxBins = hVertex->GetXaxis()->GetNbins();
310 fVertexBins.SetName("VertexBins");
311 fVertexBins.GetXaxis()->Set(fNvtxBins,-1*fZvtxCut,fZvtxCut);
312
313 TList* listOfHits = (TList*)infile.Get("Hits");
314 TList* listOfPrim = (TList*)infile.Get("Primaries");
315
316 TH1F* hEventsSelected = (TH1F*)listOfHits->FindObject("EventsSelected");
317 TH1F* hEventsAll = (TH1F*)listOfPrim->FindObject("EventsAll");
318 fListOfHits.Add(hEventsSelected);
319 fListOfPrimaries.Add(hEventsAll);
320
321 for(Int_t det =1; det<=3;det++)
322 {
323 Int_t nRings = (det==1 ? 1 : 2);
324 for(Int_t ring = 0;ring<nRings;ring++)
325 {
326 Char_t ringChar = (ring == 0 ? 'I' : 'O');
327 TH1F* allHits = (TH1F*)listOfHits->FindObject(Form("allHits_FMD%d%c",det,ringChar));
328 TH1F* doubleHits = (TH1F*)listOfHits->FindObject(Form("DoubleHits_FMD%d%c",det,ringChar));
329 fListOfHits.Add(allHits);
330 fListOfHits.Add(doubleHits);
331 for(Int_t v=0; v<fNvtxBins;v++)
332 {
333
334 TH2F* hHits = (TH2F*)listOfHits->FindObject(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v));
335 fListOfHits.Add(hHits);
336 }
337 }
338 }
339 for(Int_t iring = 0; iring<2;iring++) {
340 Char_t ringChar = (iring == 0 ? 'I' : 'O');
341 for(Int_t v=0; v<fNvtxBins;v++) {
342
343 TH2F* hPrimary = (TH2F*)listOfPrim->FindObject( Form("hPrimary_FMD_%c_vtx%d",ringChar,v));
344 fListOfPrimaries.Add(hPrimary);
345
346 }
347 }
348 GenerateCorrection();
349
350 TFile fout("backgroundFromFile.root","recreate");
351 fListOfHits.Write();
352 fListOfPrimaries.Write();
353 fListOfCorrection.Write();
354 fVertexBins.Write();
355
356 fout.Close();
357
358 if(storeInOCDB) {
359 TFile fbg("$ALICE_ROOT/FMD/Correction/Background/background.root","RECREATE");
360 fBackground->Write(AliFMDAnaParameters::GetBackgroundID());
361 fbg.Close();
362 TFile feselect("$ALICE_ROOT/FMD/Correction/EventSelectionEfficiency/eventselectionefficiency.root","RECREATE");
363 fEventSelectionEff->Write(AliFMDAnaParameters::GetEventSelectionEffID());
364 feselect.Close();
365
366 }
367
368
369
370
371}
372//_____________________________________________________________________
373//
374// EOF
375//