]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG2/FORWARD/analysis/AliFMDAnalysisTaskGenerateCorrection.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Minor fixes
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / 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"
04f1ff3d 16#include "AliMultiplicity.h"
d31dce33 17#include "AliHeader.h"
18#include "AliFMDAnaCalibBackgroundCorrection.h"
19#include "AliFMDAnaCalibEventSelectionEfficiency.h"
059c7c6b 20#include "AliGenPythiaEventHeader.h"
d31dce33 21//#include "AliCDBManager.h"
22//#include "AliCDBId.h"
23//#include "AliCDBMetaData.h"
24#include "TSystem.h"
25#include "TROOT.h"
26#include "TAxis.h"
507687cd 27#include "AliGenDPMjetEventHeader.h"
3d4a1473 28#include <TFile.h>
507687cd 29
d31dce33 30ClassImp(AliFMDAnalysisTaskGenerateCorrection)
1b418b63 31#if 0
32; // This is for emacs - do not delete
33#endif
d31dce33 34
35//_____________________________________________________________________
36AliFMDAnalysisTaskGenerateCorrection::AliFMDAnalysisTaskGenerateCorrection():
37AliAnalysisTaskSE(),
38 fListOfHits(),
39 fListOfPrimaries(),
40 fListOfCorrection(),
41 fVertexBins(),
42 fLastTrackByStrip(0),
43 fHitsByStrip(0),
44 fZvtxCut(10),
45 fNvtxBins(10),
46 fNbinsEta(200),
9f55be54 47 fBackground(0),
48 fEventSelectionEff(0)
d31dce33 49{
50 // Default constructor
51}
52//_____________________________________________________________________
53AliFMDAnalysisTaskGenerateCorrection::AliFMDAnalysisTaskGenerateCorrection(const char* name):
54 AliAnalysisTaskSE(name),
55 fListOfHits(),
56 fListOfPrimaries(),
57 fListOfCorrection(),
58 fVertexBins(),
59 fLastTrackByStrip(0),
60 fHitsByStrip(0),
61 fZvtxCut(10),
62 fNvtxBins(10),
63 fNbinsEta(200),
9f55be54 64 fBackground(0),
65 fEventSelectionEff(0)
d31dce33 66{
67
68 DefineOutput(1, TList::Class());
69 DefineOutput(2, TList::Class());
70 DefineOutput(3, TH1F::Class());
71 DefineOutput(4, TList::Class());
72}
73//_____________________________________________________________________
74void AliFMDAnalysisTaskGenerateCorrection::UserCreateOutputObjects()
75{
76// Create the output containers
77//
78
79 std::cout<<"Creating output objects"<<std::endl;
507687cd 80
df2a9c32 81 // SelectCollisionCandidates(AliVEvent::kAny);
507687cd 82
04f1ff3d 83 for(Int_t v=0; v<fNvtxBins;v++) {
84
85 TH2F* hSPDhits = new TH2F(Form("hSPDhits_vtx%d",v),
86 Form("hSPDhits_vtx%d",v),
820f4fcc 87 fNbinsEta, -4,6, 20, 0,2*TMath::Pi());
04f1ff3d 88 hSPDhits->Sumw2();
89 fListOfHits.Add(hSPDhits);
90
f7356393 91 TH1F* hReadChannels = new TH1F(Form("hFMDReadChannels_vtx%d",v),
92 Form("hFMDReadChannels_vtx%d",v),
93 fNbinsEta,-4,6);
94 hReadChannels->Sumw2();
95 fListOfHits.Add(hReadChannels);
96 TH1F* hAllChannels = new TH1F(Form("hFMDAllChannels_vtx%d",v),
97 Form("hFMDAllChannels_vtx%d",v),
98 fNbinsEta,-4,6);
99 hAllChannels->Sumw2();
100 fListOfHits.Add(hAllChannels);
101
102
103
04f1ff3d 104 for(Int_t iring = 0; iring<2;iring++) {
105 Char_t ringChar = (iring == 0 ? 'I' : 'O');
106 Int_t nSec = (iring == 1 ? 40 : 20);
107
d31dce33 108 TH2F* hPrimary = new TH2F(Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
109 Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
820f4fcc 110 fNbinsEta, -4,6, nSec, 0,2*TMath::Pi());
d31dce33 111 hPrimary->Sumw2();
112 fListOfPrimaries.Add(hPrimary);
507687cd 113 TH2F* hPrimaryNSD = new TH2F(Form("hPrimaryNSD_FMD_%c_vtx%d",ringChar,v),
114 Form("hPrimaryNSD_FMD_%c_vtx%d",ringChar,v),
115 fNbinsEta, -4,6, nSec, 0,2*TMath::Pi());
116 hPrimaryNSD->Sumw2();
117 fListOfPrimaries.Add(hPrimaryNSD);
118
04f1ff3d 119
120
d31dce33 121 }
122 }
123
04f1ff3d 124
d31dce33 125
126
127 for(Int_t det =1; det<=3;det++) {
128 Int_t nRings = (det==1 ? 1 : 2);
129 for(Int_t ring = 0;ring<nRings;ring++) {
130 Int_t nSec = (ring == 1 ? 40 : 20);
131 Char_t ringChar = (ring == 0 ? 'I' : 'O');
132 TH1F* doubleHits = new TH1F(Form("DoubleHits_FMD%d%c",det,ringChar),
820f4fcc 133 Form("DoubleHits_FMD%d%c",det,ringChar),fNbinsEta, -4,6);
d31dce33 134 TH1F* allHits = new TH1F(Form("allHits_FMD%d%c",det,ringChar),
820f4fcc 135 Form("allHits_FMD%d%c",det,ringChar), fNbinsEta, -4,6);
d31dce33 136
137 doubleHits->Sumw2();
138 allHits->Sumw2();
139 fListOfHits.Add(allHits);
140 fListOfHits.Add(doubleHits);
141
142 for(Int_t v=0; v<fNvtxBins;v++) {
143 TH2F* hHits = new TH2F(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
144 Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
820f4fcc 145 fNbinsEta, -4,6, nSec, 0, 2*TMath::Pi());
d31dce33 146 hHits->Sumw2();
147 fListOfHits.Add(hHits);
507687cd 148 TH2F* hHitsNSD = new TH2F(Form("hHitsNSD_FMD%d%c_vtx%d", det,ringChar,v),
149 Form("hHitsNSD_FMD%d%c_vtx%d", det,ringChar,v),
150 fNbinsEta, -4,6, nSec, 0, 2*TMath::Pi());
151 hHitsNSD->Sumw2();
152 fListOfHits.Add(hHitsNSD);
153
d31dce33 154
155 }
156 }
157 }
158
159 TH1F* hEventsSelected = new TH1F("EventsSelected","EventsSelected",fNvtxBins,0,fNvtxBins);
160 TH1F* hEventsAll = new TH1F("EventsAll","EventsAll",fNvtxBins,0,fNvtxBins);
059c7c6b 161 TH1F* hEventsAllNSD = new TH1F("EventsAllNSD","EventsAllNSD",fNvtxBins,0,fNvtxBins);
f55d559b 162 TH1F* hEventsSelectedVtx = new TH1F("EventsSelectedVtx","EventsSelectedVtx",fNvtxBins,0,fNvtxBins);
163 TH1F* hEventsSelectedTrigger = new TH1F("EventsSelectedTrigger","EventsSelectedTrigger",fNvtxBins,0,fNvtxBins);
f55d559b 164
059c7c6b 165 TH1F* hEventsSelectedNSDVtx = new TH1F("EventsSelectedNSDVtx","EventsSelectedNSDVtx",fNvtxBins,0,fNvtxBins);
166 TH1F* hEventsSelectedNSD = new TH1F("EventsSelectedNSD","EventsSelectedNSD",fNvtxBins,0,fNvtxBins);
167
850ad94d 168 TH1F* hXvtx = new TH1F("hXvtx","x vertex distribution",100,-2,2);
169 TH1F* hYvtx = new TH1F("hYvtx","y vertex distribution",100,-2,2);
170 TH1F* hZvtx = new TH1F("hZvtx","z vertex distribution",4*fNvtxBins,-4*fZvtxCut,4*fZvtxCut);
171
172 fListOfPrimaries.Add(hXvtx);
173 fListOfPrimaries.Add(hYvtx);
174 fListOfPrimaries.Add(hZvtx);
f55d559b 175
0b0a4ae5 176 hEventsSelected->Sumw2();
177 hEventsAll->Sumw2();
059c7c6b 178 hEventsAllNSD->Sumw2();
179
d31dce33 180 fListOfHits.Add(hEventsSelected);
f55d559b 181 fListOfHits.Add(hEventsSelectedVtx);
182 fListOfHits.Add(hEventsSelectedTrigger);
059c7c6b 183 fListOfHits.Add(hEventsSelectedNSDVtx);
184 fListOfHits.Add(hEventsSelectedNSD);
185
186
d31dce33 187 fListOfPrimaries.Add(hEventsAll);
059c7c6b 188 fListOfPrimaries.Add(hEventsAllNSD);
d31dce33 189
ab3e0abc 190 for(Int_t v=0; v<fNvtxBins;v++) {
191
192 for(Int_t iring = 0; iring<2;iring++) {
193 Char_t ringChar = (iring == 0 ? 'I' : 'O');
194 Int_t nSec = (iring == 1 ? 40 : 20);
195 TH2F* hAnalysed = new TH2F(Form("Analysed_FMD%c_vtx%d",ringChar,v),
196 Form("Analysed_FMD%c_vtx%d",ringChar,v),
820f4fcc 197 fNbinsEta, -4,6, nSec, 0,2*TMath::Pi());
ab3e0abc 198
199 hAnalysed->Sumw2();
200 fListOfPrimaries.Add(hAnalysed);
201
059c7c6b 202 TH2F* hAnalysedNSD = new TH2F(Form("AnalysedNSD_FMD%c_vtx%d",ringChar,v),
203 Form("AnalysedNSD_FMD%c_vtx%d",ringChar,v),
204 fNbinsEta, -4,6, nSec, 0,2*TMath::Pi());
205
206 hAnalysedNSD->Sumw2();
207 fListOfPrimaries.Add(hAnalysedNSD);
208
ab3e0abc 209 TH2F* hInel = new TH2F(Form("Inel_FMD%c_vtx%d",ringChar,v),
210 Form("Inel_FMD%c_vtx%d",ringChar,v),
820f4fcc 211 fNbinsEta, -4,6, nSec, 0,2*TMath::Pi());
ab3e0abc 212
213 hInel->Sumw2();
214 fListOfPrimaries.Add(hInel);
215
059c7c6b 216 TH2F* hNSD = new TH2F(Form("NSD_FMD%c_vtx%d",ringChar,v),
217 Form("NSD_FMD%c_vtx%d",ringChar,v),
218 fNbinsEta, -4,6, nSec, 0,2*TMath::Pi());
219
220 hNSD->Sumw2();
221 fListOfPrimaries.Add(hNSD);
222
ab3e0abc 223 }
224 }
d31dce33 225
226 fVertexBins.SetName("VertexBins");
227 fVertexBins.GetXaxis()->Set(fNvtxBins,-1*fZvtxCut,fZvtxCut);
228
507687cd 229
230
d31dce33 231}
232//_____________________________________________________________________
233void AliFMDAnalysisTaskGenerateCorrection::Init()
234{
235 fLastTrackByStrip.Reset(-1);
236
237
238}
239//_____________________________________________________________________
240void AliFMDAnalysisTaskGenerateCorrection::UserExec(Option_t */*option*/)
241{
242
243 fLastTrackByStrip.Reset(-1);
244 fHitsByStrip.Reset(0);
245 AliMCEvent* mcevent = MCEvent();
246
247 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
248
249
250 AliESDEvent* esdevent = (AliESDEvent*)InputEvent();
059c7c6b 251
252 pars->SetTriggerStatus(esdevent);
253
d31dce33 254 Double_t esdvertex[3];
f55d559b 255 Bool_t vtxStatus = pars->GetVertex(esdevent,esdvertex);
d31dce33 256
f7356393 257 /* Double_t deltaEsd = 2*fZvtxCut/fNvtxBins;
258 Double_t vertexBinDoubleEsd = (esdvertex[2] + fZvtxCut) / deltaEsd;
259 Int_t vertexBinEsd = (Int_t)vertexBinDoubleEsd;*/
260
d31dce33 261 AliMCParticle* particle = 0;
262 AliStack* stack = mcevent->Stack();
263
264 UShort_t det,sec,strip;
265 Char_t ring;
266
267 Int_t nTracks = mcevent->GetNumberOfTracks();
268 AliHeader* header = mcevent->Header();
269 AliGenEventHeader* genHeader = header->GenEventHeader();
270
059c7c6b 271 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
507687cd 272 AliGenDPMjetEventHeader* dpmHeader = dynamic_cast<AliGenDPMjetEventHeader*>(header->GenEventHeader());
273
274 if (!pythiaGenHeader && !dpmHeader) {
275 std::cout<<" no pythia or dpm header! - NSD selection unusable"<<std::endl;
04d7f248 276 //return;
059c7c6b 277 }
059c7c6b 278 Bool_t nsd = kTRUE;
04d7f248 279 if(pythiaGenHeader) {
280 Int_t pythiaType = pythiaGenHeader->ProcessType();
281
282 if(pythiaType==92 || pythiaType==93)
283 nsd = kFALSE;
284 }
507687cd 285 if(dpmHeader) {
286 Int_t processType = dpmHeader->ProcessType();
287 if(processType == 5 || processType == 6)
288 nsd = kFALSE;
289
290 }
291
059c7c6b 292
d31dce33 293 TArrayF vertex;
294 genHeader->PrimaryVertex(vertex);
295
850ad94d 296 TH1F* hXvtx = (TH1F*)fListOfPrimaries.FindObject("hXvtx");
297 hXvtx->Fill(vertex.At(0));
298 TH1F* hYvtx = (TH1F*)fListOfPrimaries.FindObject("hYvtx");
299 hYvtx->Fill(vertex.At(1));
300 TH1F* hZvtx = (TH1F*)fListOfPrimaries.FindObject("hZvtx");
301 hZvtx->Fill(vertex.At(2));
302
303
059c7c6b 304
d31dce33 305
306 Double_t delta = 2*fZvtxCut/fNvtxBins;
307 Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
308 Int_t vertexBin = (Int_t)vertexBinDouble;
309
310 // Vertex determination correction
f55d559b 311 TH1F* hEventsSelected = (TH1F*)fListOfHits.FindObject("EventsSelected");
312 TH1F* hEventsSelectedVtx = (TH1F*)fListOfHits.FindObject("EventsSelectedVtx");
313 TH1F* hEventsSelectedTrigger = (TH1F*)fListOfHits.FindObject("EventsSelectedTrigger");
059c7c6b 314 TH1F* hEventsSelectedNSDVtx = (TH1F*)fListOfHits.FindObject("EventsSelectedNSDVtx");
315 TH1F* hEventsSelectedNSD = (TH1F*)fListOfHits.FindObject("EventsSelectedNSD");
f55d559b 316 TH1F* hEventsAll = (TH1F*)fListOfPrimaries.FindObject("EventsAll");
059c7c6b 317 TH1F* hEventsAllNSD = (TH1F*)fListOfPrimaries.FindObject("EventsAllNSD");
f55d559b 318
d31dce33 319 // TH1F* hTriggered = (TH1F*)fListOfHits.FindObject("Triggered");
320 // TH1F* hTriggeredAll = (TH1F*)fListOfPrimaries.FindObject("TriggeredAll");
321
322 Bool_t vtxFound = kTRUE;
f55d559b 323 if(!vtxStatus)
d31dce33 324 vtxFound = kFALSE;
d31dce33 325
4b8bdb60 326 //if(TMath::Abs(vertex.At(2)) > fZvtxCut) {
327 // vtxFound = kFALSE;
059c7c6b 328
4b8bdb60 329 //}
df2a9c32 330
331 if(TMath::Abs(vertex.At(2)) > fZvtxCut) {
332 return;
333 }
334
335
059c7c6b 336 Bool_t isTriggered = pars->IsEventTriggered(AliFMDAnaParameters::kMB1);
337 Bool_t isTriggeredNSD = pars->IsEventTriggered(AliFMDAnaParameters::kNSD);
f55d559b 338 if(vtxFound && isTriggered) hEventsSelected->Fill(vertexBin);
d31dce33 339
f55d559b 340 if(vtxFound) hEventsSelectedVtx->Fill(vertexBin);
341 if(isTriggered) hEventsSelectedTrigger->Fill(vertexBin);
059c7c6b 342
343 if(vtxFound && isTriggeredNSD) hEventsSelectedNSDVtx->Fill(vertexBin);
344 if(isTriggeredNSD) hEventsSelectedNSD->Fill(vertexBin);
345
346
d31dce33 347 hEventsAll->Fill(vertexBin);
059c7c6b 348 if(nsd) hEventsAllNSD->Fill(vertexBin);
70d74659 349
f7356393 350 // FMD dead channels
351
352 TH1F* hFMDReadChannels = (TH1F*)fListOfHits.FindObject(Form("hFMDReadChannels_vtx%d",vertexBin));
353 TH1F* hFMDAllChannels = (TH1F*)fListOfHits.FindObject(Form("hFMDAllChannels_vtx%d",vertexBin));
354
355 AliESDFMD* fmd = esdevent->GetFMDData();
356
357 for(UShort_t d=1;d<=3;d++) {
358 Int_t nRings = (d==1 ? 1 : 2);
359 for (UShort_t ir = 0; ir < nRings; ir++) {
360 Char_t r = (ir == 0 ? 'I' : 'O');
361 UShort_t nsec = (ir == 0 ? 20 : 40);
362 UShort_t nstr = (ir == 0 ? 512 : 256);
363 for(UShort_t s =0; s < nsec; s++) {
364 for(UShort_t str = 0; str < nstr; str++) {
365 Float_t mult = fmd->Multiplicity(d,r,s,str);
366
367 hFMDAllChannels->Fill(pars->GetEtaFromStrip(d,r,s,str,esdvertex[2]));
368
369 if(mult != AliESDFMD::kInvalidMult)
370 hFMDReadChannels->Fill(pars->GetEtaFromStrip(d,r,s,str,esdvertex[2]));
371 }
372 }
373 }
374 }
375 // End, FMD dead channels
df2a9c32 376
d31dce33 377 for(Int_t i = 0 ;i<nTracks;i++) {
7aad0c47 378 particle = (AliMCParticle*) mcevent->GetTrack(i);
d31dce33 379
380 if(!particle)
381 continue;
382
383 if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
384
df2a9c32 385 if(vtxFound && isTriggered) {
386 TH2F* hPrimaryInner = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",'I',vertexBin));
387 TH2F* hPrimaryOuter = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",'O',vertexBin));
388 hPrimaryInner->Fill(particle->Eta(),particle->Phi());
389 hPrimaryOuter->Fill(particle->Eta(),particle->Phi());
390 if( isTriggeredNSD) {
391 TH2F* hPrimaryInnerNSD = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimaryNSD_FMD_%c_vtx%d",'I',vertexBin));
392 TH2F* hPrimaryOuterNSD = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimaryNSD_FMD_%c_vtx%d",'O',vertexBin));
393 hPrimaryInnerNSD->Fill(particle->Eta(),particle->Phi());
394 hPrimaryOuterNSD->Fill(particle->Eta(),particle->Phi());
395 }
507687cd 396 }
ab3e0abc 397 TH2F* hAnalysedInner = (TH2F*)fListOfPrimaries.FindObject( Form("Analysed_FMD%c_vtx%d",'I',vertexBin));
398 TH2F* hAnalysedOuter = (TH2F*)fListOfPrimaries.FindObject( Form("Analysed_FMD%c_vtx%d",'O',vertexBin));
059c7c6b 399 TH2F* hAnalysedNSDInner = (TH2F*)fListOfPrimaries.FindObject( Form("AnalysedNSD_FMD%c_vtx%d",'I',vertexBin));
400 TH2F* hAnalysedNSDOuter = (TH2F*)fListOfPrimaries.FindObject( Form("AnalysedNSD_FMD%c_vtx%d",'O',vertexBin));
ab3e0abc 401 TH2F* hInelInner = (TH2F*)fListOfPrimaries.FindObject( Form("Inel_FMD%c_vtx%d",'I',vertexBin));
402 TH2F* hInelOuter = (TH2F*)fListOfPrimaries.FindObject( Form("Inel_FMD%c_vtx%d",'O',vertexBin));
059c7c6b 403 TH2F* hNSDInner = (TH2F*)fListOfPrimaries.FindObject( Form("NSD_FMD%c_vtx%d",'I',vertexBin));
404 TH2F* hNSDOuter = (TH2F*)fListOfPrimaries.FindObject( Form("NSD_FMD%c_vtx%d",'O',vertexBin));
ab3e0abc 405
820f4fcc 406 //if(vtxFound && isTriggered) {
059c7c6b 407 if(isTriggeredNSD) {
408 hAnalysedNSDInner->Fill(particle->Eta(),particle->Phi());
409 hAnalysedNSDOuter->Fill(particle->Eta(),particle->Phi());
410 }
df2a9c32 411 if(isTriggered && vtxFound) {
ab3e0abc 412 hAnalysedInner->Fill(particle->Eta(),particle->Phi());
413 hAnalysedOuter->Fill(particle->Eta(),particle->Phi());
414 }
415 hInelInner->Fill(particle->Eta(),particle->Phi());
416 hInelOuter->Fill(particle->Eta(),particle->Phi());
417
059c7c6b 418 if(nsd) {
419 hNSDInner->Fill(particle->Eta(),particle->Phi());
420 hNSDOuter->Fill(particle->Eta(),particle->Phi());
421 }
d31dce33 422 }
df2a9c32 423 if(!vtxFound || !isTriggered) continue;
d31dce33 424
425 for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) {
426
427 AliTrackReference* ref = particle->GetTrackReference(j);
428
429 if(ref->DetectorId() != AliTrackReference::kFMD)
430 continue;
ab3e0abc 431
d31dce33 432 AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip);
433 Float_t thisStripTrack = fLastTrackByStrip(det,ring,sec,strip);
434 if(particle->Charge() != 0 && i != thisStripTrack ) {
435
436 Float_t phi = pars->GetPhiFromSector(det,ring,sec);
437 Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));
438
439 TH2F* hHits = (TH2F*)fListOfHits.FindObject(Form("hHits_FMD%d%c_vtx%d", det,ring,vertexBin));
440 hHits->Fill(eta,phi);
507687cd 441 if( isTriggeredNSD) {
442 TH2F* hHitsNSD = (TH2F*)fListOfHits.FindObject(Form("hHitsNSD_FMD%d%c_vtx%d", det,ring,vertexBin));
443 hHitsNSD->Fill(eta,phi);
444 }
445
d31dce33 446 Float_t nstrips = (ring =='O' ? 256 : 512);
447 fHitsByStrip(det,ring,sec,strip) +=1;
448 TH1F* allHits = (TH1F*)fListOfHits.FindObject(Form("allHits_FMD%d%c",det,ring));
449 TH1F* doubleHits = (TH1F*)fListOfHits.FindObject(Form("DoubleHits_FMD%d%c",det,ring));
450
451 if(fHitsByStrip(det,ring,sec,strip) == 1)
452 allHits->Fill(eta);
453
454 doubleHits->Fill(eta);
455
456 fLastTrackByStrip(det,ring,sec,strip) = (Float_t)i;
457 if(strip >0)
458 fLastTrackByStrip(det,ring,sec,strip-1) = (Float_t)i;
459 if(strip < (nstrips - 1))
460 fLastTrackByStrip(det,ring,sec,strip+1) = (Float_t)i;
461 }
462 }
463
464 }
04f1ff3d 465
04f1ff3d 466
df2a9c32 467 if(vtxFound && isTriggered) {
468 //SPD part HHD
469 TH2F* hSPDMult = (TH2F*)fListOfHits.FindObject(Form("hSPDhits_vtx%d", vertexBin));
470
471 const AliMultiplicity* spdmult = esdevent->GetMultiplicity();
472 for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++)
473 hSPDMult->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
d31dce33 474
df2a9c32 475 for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++)
476 hSPDMult->Fill(-TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.)),spdmult->GetPhiSingle(j));
477
478 }
479
d31dce33 480 PostData(1, &fListOfHits);
481 PostData(2, &fListOfPrimaries);
482 PostData(3, &fVertexBins);
483}
484//_____________________________________________________________________
485void AliFMDAnalysisTaskGenerateCorrection::Terminate(Option_t */*option*/)
486{
487 /* TH1F* allHits = (TH1F*)fListOfHits.FindObject("allHits");
488 TH1F* doubleHits = (TH1F*)fListOfHits.FindObject("DoubleHits");
489
490 doubleHits->Divide(allHits);
491 GenerateCorrection();
492 PostData(1, &fListOfHits);
493 PostData(4, &fListOfCorrection);*/
494
495}
496//_____________________________________________________________________
497void AliFMDAnalysisTaskGenerateCorrection::GenerateCorrection() {
498
499 fBackground = new AliFMDAnaCalibBackgroundCorrection();
500 fEventSelectionEff = new AliFMDAnaCalibEventSelectionEfficiency();
501
ab3e0abc 502 //Event selection
f55d559b 503 TH1F* hEventsSelected = (TH1F*)fListOfHits.FindObject("EventsSelected");
504 TH1F* hEventsSelectedVtx = (TH1F*)fListOfHits.FindObject("EventsSelectedVtx");
505 TH1F* hEventsSelectedTrigger = (TH1F*)fListOfHits.FindObject("EventsSelectedTrigger");
059c7c6b 506 TH1F* hEventsSelectedNSDVtx = (TH1F*)fListOfHits.FindObject("EventsSelectedNSDVtx");
507 TH1F* hEventsSelectedNSD = (TH1F*)fListOfHits.FindObject("EventsSelectedNSD");
508
f55d559b 509 TH1F* hEventsAll = (TH1F*)fListOfPrimaries.FindObject("EventsAll");
059c7c6b 510 TH1F* hEventsAllNSD = (TH1F*)fListOfPrimaries.FindObject("EventsAllNSD");
820f4fcc 511 TH1F* hEventsSelectedVtxDivByTr = (TH1F*)hEventsSelectedVtx->Clone("hEventsSelectedVtxDivByTr");
059c7c6b 512 TH1F* hEventsSelectedNSDVtxDivByNSD = (TH1F*)hEventsSelectedNSDVtx->Clone("hEventsSelectedNSDVtxDivByNSD");
513
820f4fcc 514 fListOfHits.Add(hEventsSelectedVtxDivByTr);
059c7c6b 515 fListOfHits.Add(hEventsSelectedNSDVtxDivByNSD);
516 hEventsSelectedNSDVtxDivByNSD->Divide(hEventsSelectedNSD);
820f4fcc 517 hEventsSelectedVtxDivByTr->Divide(hEventsSelectedTrigger);
059c7c6b 518
519 for(Int_t v=0;v<fNvtxBins ; v++) {
ab3e0abc 520 TH2F* hAnalysedInner = (TH2F*)fListOfPrimaries.FindObject(Form("Analysed_FMD%c_vtx%d",'I',v));
521 TH2F* hAnalysedOuter = (TH2F*)fListOfPrimaries.FindObject(Form("Analysed_FMD%c_vtx%d",'O',v));
059c7c6b 522 TH2F* hAnalysedNSDInner = (TH2F*)fListOfPrimaries.FindObject( Form("AnalysedNSD_FMD%c_vtx%d",'I',v));
523 TH2F* hAnalysedNSDOuter = (TH2F*)fListOfPrimaries.FindObject( Form("AnalysedNSD_FMD%c_vtx%d",'O',v));
524
ab3e0abc 525 TH2F* hInelInner = (TH2F*)fListOfPrimaries.FindObject(Form("Inel_FMD%c_vtx%d",'I',v));
526 TH2F* hInelOuter = (TH2F*)fListOfPrimaries.FindObject(Form("Inel_FMD%c_vtx%d",'O',v));
059c7c6b 527 TH2F* hNSDInner = (TH2F*)fListOfPrimaries.FindObject( Form("NSD_FMD%c_vtx%d",'I',v));
528 TH2F* hNSDOuter = (TH2F*)fListOfPrimaries.FindObject( Form("NSD_FMD%c_vtx%d",'O',v));
820f4fcc 529 // hAnalysedInner->Scale((hEventsSelected->GetBinContent(v+1) == 0 ? 0 : 1/hEventsSelected->GetBinContent(v+1)));
530 //hAnalysedOuter->Scale((hEventsSelected->GetBinContent(v+1) == 0 ? 0 : 1/hEventsSelected->GetBinContent(v+1)));
531
532 hAnalysedInner->Scale((hEventsSelectedTrigger->GetBinContent(v+1) == 0 ? 0 : 1/hEventsSelectedTrigger->GetBinContent(v+1)));
533
534 hAnalysedOuter->Scale((hEventsSelectedTrigger->GetBinContent(v+1) == 0 ? 0 : 1/hEventsSelectedTrigger->GetBinContent(v+1)));
059c7c6b 535 hAnalysedNSDInner->Scale((hEventsSelectedNSD->GetBinContent(v+1) == 0 ? 0 : 1/hEventsSelectedNSD->GetBinContent(v+1)));
536
537 hAnalysedNSDOuter->Scale((hEventsSelectedNSD->GetBinContent(v+1) == 0 ? 0 : 1/hEventsSelectedNSD->GetBinContent(v+1)));
538
820f4fcc 539
059c7c6b 540
ab3e0abc 541 hInelInner->Scale((hEventsAll->GetBinContent(v+1) == 0 ? 0 : 1/hEventsAll->GetBinContent(v+1)));
542 hInelOuter->Scale((hEventsAll->GetBinContent(v+1) == 0 ? 0 : 1/hEventsAll->GetBinContent(v+1)));
820f4fcc 543
059c7c6b 544 hNSDInner->Scale((hEventsAllNSD->GetBinContent(v+1) == 0 ? 0 : 1/hEventsAllNSD->GetBinContent(v+1)));
545 hNSDOuter->Scale((hEventsAllNSD->GetBinContent(v+1) == 0 ? 0 : 1/hEventsAllNSD->GetBinContent(v+1)));
546
820f4fcc 547
ab3e0abc 548 hAnalysedInner->Divide(hInelInner);
549 hAnalysedOuter->Divide(hInelOuter);
820f4fcc 550
059c7c6b 551 hAnalysedNSDInner->Divide(hNSDInner);
552 hAnalysedNSDOuter->Divide(hNSDOuter);
553
554 fEventSelectionEff->SetCorrection("INEL",v,'I',hAnalysedInner);
555 fEventSelectionEff->SetCorrection("INEL",v,'O',hAnalysedOuter);
556 //NSD
557 fEventSelectionEff->SetCorrection("NSD",v,'I',hAnalysedNSDInner);
558 fEventSelectionEff->SetCorrection("NSD",v,'O',hAnalysedNSDOuter);
559
f7356393 560 //FMD dead channels
561 TH1F* hFMDReadChannels = (TH1F*)fListOfHits.FindObject(Form("hFMDReadChannels_vtx%d",v));
562
563 TH1F* hFMDAllChannels = (TH1F*)fListOfHits.FindObject(Form("hFMDAllChannels_vtx%d",v));
564 if(hFMDReadChannels) {
565 TH1F* hFMDDeadCorrection = (TH1F*)hFMDReadChannels->Clone("hFMDDeadCorrection");
566 hFMDDeadCorrection->Divide(hFMDAllChannels);
567 fBackground->SetFMDDeadCorrection(v,hFMDDeadCorrection);
568 fListOfCorrection.Add(hFMDDeadCorrection);
569 }
570 else AliWarning("No Dead Channel Correction generated");
059c7c6b 571
541c19ed 572 }
ab3e0abc 573
820f4fcc 574 Float_t vtxEff = 1;
575 if(hEventsSelectedTrigger->GetEntries())
576 vtxEff = hEventsSelectedVtx->GetEntries() / hEventsSelectedTrigger->GetEntries();
577 fEventSelectionEff->SetVtxToTriggerRatio(vtxEff);
ab3e0abc 578
b85ea106 579 // hEventsAll->Divide(hEventsAll,hEventsSelected,1,1,"B");
f55d559b 580 hEventsSelectedVtx->Divide(hEventsAll);
581 hEventsSelectedTrigger->Divide(hEventsAll);
0b0a4ae5 582
059c7c6b 583 hEventsSelectedNSDVtx->Divide(hEventsAllNSD);
584 hEventsSelectedNSD->Divide(hEventsAllNSD);
585
0b0a4ae5 586 for(Int_t i = 1; i<=hEventsSelected->GetNbinsX(); i++) {
587 if(hEventsSelected->GetBinContent(i) == 0 )
588 continue;
059c7c6b 589 Float_t b = hEventsSelected->GetBinContent(i);
590 // Float_t b = hEventsSelected->GetBinContent(i)*hEventsSelectedTrigger->GetBinContent(i);
b85ea106 591 Float_t db = hEventsSelected->GetBinError(i);
0b0a4ae5 592 Float_t sum = hEventsAll->GetBinContent(i);
593 Float_t dsum = hEventsAll->GetBinError(i);
b85ea106 594 Float_t a = sum-b;
595 Float_t da = TMath::Sqrt(TMath::Power(db,2) + TMath::Power(dsum,2));
0b0a4ae5 596
b85ea106 597 Float_t cor = sum / b;
598 Float_t ecor = TMath::Sqrt(TMath::Power(da,2) + TMath::Power(a/(b*db),2)) / b;
0b0a4ae5 599
b85ea106 600 hEventsAll->SetBinContent(i,cor);
601 hEventsAll->SetBinError(i,ecor);
0b0a4ae5 602
603 }
d31dce33 604
ab3e0abc 605 //TH1F* hEventTest = (TH1F*)hEventsAll->Clone("hEventTest");
606
607
608
b85ea106 609 fEventSelectionEff->SetCorrection(hEventsAll);
d31dce33 610
611 for(Int_t det= 1; det <=3; det++) {
612 Int_t nRings = (det==1 ? 1 : 2);
613
614 for(Int_t iring = 0; iring<nRings; iring++) {
615 Char_t ring = (iring == 0 ? 'I' : 'O');
616 TH1F* allHits = (TH1F*)fListOfHits.FindObject(Form("allHits_FMD%d%c",det,ring));
617 TH1F* doubleHits = (TH1F*)fListOfHits.FindObject(Form("DoubleHits_FMD%d%c",det,ring));
b85ea106 618 allHits->Divide(doubleHits);
619
620 fBackground->SetDoubleHitCorrection(det,ring,allHits);
621
d31dce33 622 for(Int_t vertexBin=0;vertexBin<fNvtxBins ;vertexBin++) {
623 TH2F* hHits = (TH2F*)fListOfHits.FindObject(Form("hHits_FMD%d%c_vtx%d", det,ring,vertexBin));
507687cd 624 TH2F* hPrimary = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",ring,vertexBin));
625 TH2F* hCorrection = (TH2F*)hHits->Clone(Form("FMD%d%c_vtxbin_%d_correction",det,ring,vertexBin));
d31dce33 626 hCorrection->Divide(hPrimary);
507687cd 627
628 TH2F* hHitsNSD = (TH2F*)fListOfHits.FindObject(Form("hHitsNSD_FMD%d%c_vtx%d", det,ring,vertexBin));
629 TH2F* hPrimaryNSD = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimaryNSD_FMD_%c_vtx%d",ring,vertexBin));
df2a9c32 630 TH2F* hCorrectionNSD = 0;
631 if(hHitsNSD && hPrimaryNSD) {
632 hCorrectionNSD = (TH2F*)hHitsNSD->Clone(Form("FMDNSD%d%c_vtxbin_%d_correction",det,ring,vertexBin));
633 hCorrectionNSD->Divide(hPrimaryNSD);
634 }
0b0a4ae5 635
d31dce33 636 hCorrection->SetTitle(hCorrection->GetName());
637 fListOfCorrection.Add(hCorrection);
638 fBackground->SetBgCorrection(det,ring,vertexBin,hCorrection);
df2a9c32 639
640 if(hCorrectionNSD) {
641 fBackground->SetNSDBgCorrection(det,ring,vertexBin,hCorrectionNSD);
642 hCorrectionNSD->SetTitle(hCorrectionNSD->GetName());
643 fListOfCorrection.Add(hCorrectionNSD);
644 }
d31dce33 645
646 }
647
648 }
649 }
04f1ff3d 650 for(Int_t vertexBin=0;vertexBin<fNvtxBins ;vertexBin++) {
651 TH2F* hPrimary = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",'I',vertexBin));
652 TH2F* hSPDMult = (TH2F*)fListOfHits.FindObject(Form("hSPDhits_vtx%d", vertexBin));
11e702b3 653 if(!hSPDMult) continue;
654
04f1ff3d 655 TH2F* hCorrection = (TH2F*)hSPDMult->Clone(Form("SPD_vtxbin_%d_correction",vertexBin));
656 hCorrection->SetTitle(hCorrection->GetName());
657 fListOfCorrection.Add(hCorrection);
658 hCorrection->Divide(hPrimary);
f7356393 659
04f1ff3d 660
661 TH1F* hAlive = new TH1F(Form("hAliveSPD_vtxbin%d",vertexBin),Form("hAliveSPD_vtxbin%d",vertexBin),hSPDMult->GetNbinsX(),hSPDMult->GetXaxis()->GetXmin(), hSPDMult->GetXaxis()->GetXmax());
662 TH1F* hPresent = new TH1F(Form("hPresentSPD_vtxbin%d",vertexBin),Form("hPresentSPD_vtxbin%d",vertexBin),hSPDMult->GetNbinsX(),hSPDMult->GetXaxis()->GetXmin(), hSPDMult->GetXaxis()->GetXmax());
663 for(Int_t xx = 1; xx <=hSPDMult->GetNbinsX(); xx++) {
664
df2a9c32 665
666
04f1ff3d 667 for(Int_t yy = 1; yy <=hSPDMult->GetNbinsY(); yy++) {
df2a9c32 668 if(TMath::Abs(hCorrection->GetXaxis()->GetBinCenter(xx)) > 1.9) {
669 hCorrection->SetBinContent(xx,yy,0.);
670 hCorrection->SetBinError(xx,yy,0.);
671 }
672
673 if(hCorrection->GetBinContent(xx,yy) > 0.9) {
674 hAlive->Fill(hCorrection->GetXaxis()->GetBinCenter(xx) );
675 }
676 else {
677 hCorrection->SetBinContent(xx,yy,0.);
678 hCorrection->SetBinError(xx,yy,0.);
679 }
680
04f1ff3d 681 hPresent->Fill(hCorrection->GetXaxis()->GetBinCenter(xx));
682
683 }
684 }
685 TH1F* hDeadCorrection = (TH1F*)hAlive->Clone(Form("hSPDDeadCorrection_vtxbin%d",vertexBin));
686 hDeadCorrection->Divide(hPresent);
687 fBackground->SetSPDDeadCorrection(vertexBin,hDeadCorrection);
688 fListOfCorrection.Add(hDeadCorrection);
ab3e0abc 689
f7356393 690 fBackground->SetBgCorrection(0,'Q',vertexBin,hCorrection);
691 }
ab3e0abc 692
693
d31dce33 694 TAxis refAxis(fNvtxBins,-1*fZvtxCut,fZvtxCut);
695 fBackground->SetRefAxis(&refAxis);
696
697}
698//_____________________________________________________________________
699void AliFMDAnalysisTaskGenerateCorrection::ReadFromFile(const Char_t* filename, Bool_t storeInOCDB, Int_t /*runNo*/) {
0b0a4ae5 700
d31dce33 701 TFile infile(filename);
702 TH1F* hVertex = (TH1F*)infile.Get("VertexBins");
703 fZvtxCut = hVertex->GetXaxis()->GetXmax();
704 fNvtxBins = hVertex->GetXaxis()->GetNbins();
705 fVertexBins.SetName("VertexBins");
706 fVertexBins.GetXaxis()->Set(fNvtxBins,-1*fZvtxCut,fZvtxCut);
707
708 TList* listOfHits = (TList*)infile.Get("Hits");
709 TList* listOfPrim = (TList*)infile.Get("Primaries");
710
f55d559b 711 TH1F* hEventsSelected = (TH1F*)listOfHits->FindObject("EventsSelected");
712 TH1F* hEventsSelectedVtx = (TH1F*)listOfHits->FindObject("EventsSelectedVtx");
713 TH1F* hEventsSelectedTrigger = (TH1F*)listOfHits->FindObject("EventsSelectedTrigger");
059c7c6b 714 TH1F* hEventsSelectedNSDVtx = (TH1F*)listOfHits->FindObject("EventsSelectedNSDVtx");
715 TH1F* hEventsSelectedNSD = (TH1F*)listOfHits->FindObject("EventsSelectedNSD");
f55d559b 716 TH1F* hEventsAll = (TH1F*)listOfPrim->FindObject("EventsAll");
059c7c6b 717 TH1F* hEventsAllNSD = (TH1F*)listOfPrim->FindObject("EventsAllNSD");
f55d559b 718
d31dce33 719 fListOfHits.Add(hEventsSelected);
f55d559b 720 fListOfHits.Add(hEventsSelectedVtx);
059c7c6b 721 fListOfHits.Add(hEventsSelectedNSD);
722 fListOfHits.Add(hEventsSelectedNSDVtx);
f55d559b 723 fListOfHits.Add(hEventsSelectedTrigger);
d31dce33 724 fListOfPrimaries.Add(hEventsAll);
059c7c6b 725 fListOfPrimaries.Add(hEventsAllNSD);
d31dce33 726
850ad94d 727 TH1F* hXvtx = (TH1F*)listOfPrim->FindObject("hXvtx");
728 TH1F* hYvtx = (TH1F*)listOfPrim->FindObject("hYvtx");
729 TH1F* hZvtx = (TH1F*)listOfPrim->FindObject("hZvtx");
730 fListOfPrimaries.Add(hXvtx);
731 fListOfPrimaries.Add(hYvtx);
732 fListOfPrimaries.Add(hZvtx);
733
d31dce33 734 for(Int_t det =1; det<=3;det++)
059c7c6b 735 {
736 Int_t nRings = (det==1 ? 1 : 2);
737 for(Int_t ring = 0;ring<nRings;ring++)
738 {
739 Char_t ringChar = (ring == 0 ? 'I' : 'O');
740 TH1F* allHits = (TH1F*)listOfHits->FindObject(Form("allHits_FMD%d%c",det,ringChar));
741 TH1F* doubleHits = (TH1F*)listOfHits->FindObject(Form("DoubleHits_FMD%d%c",det,ringChar));
742 fListOfHits.Add(allHits);
743 fListOfHits.Add(doubleHits);
744 for(Int_t v=0; v<fNvtxBins;v++)
745 {
746
747 TH2F* hHits = (TH2F*)listOfHits->FindObject(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v));
748 fListOfHits.Add(hHits);
507687cd 749 TH2F* hHitsNSD = (TH2F*)listOfHits->FindObject(Form("hHitsNSD_FMD%d%c_vtx%d", det,ringChar,v));
df2a9c32 750 if(hHitsNSD)
751 fListOfHits.Add(hHitsNSD);
059c7c6b 752 }
753 }
754 }
04f1ff3d 755 for(Int_t v=0; v<fNvtxBins;v++) {
756 TH2F* hSPDHits = (TH2F*)listOfHits->FindObject(Form("hSPDhits_vtx%d", v));
757 fListOfHits.Add(hSPDHits);
f7356393 758 TH1F* hFMDReadChannels = (TH1F*)listOfHits->FindObject(Form("hFMDReadChannels_vtx%d",v));
759 TH1F* hFMDAllChannels = (TH1F*)listOfHits->FindObject(Form("hFMDAllChannels_vtx%d",v));
04f1ff3d 760
f7356393 761 if(hFMDReadChannels && hFMDAllChannels) {
762 fListOfHits.Add(hFMDReadChannels);
763 fListOfHits.Add(hFMDAllChannels);
764 }
04f1ff3d 765 for(Int_t iring = 0; iring<2;iring++) {
766 Char_t ringChar = (iring == 0 ? 'I' : 'O');
767
d31dce33 768
769 TH2F* hPrimary = (TH2F*)listOfPrim->FindObject( Form("hPrimary_FMD_%c_vtx%d",ringChar,v));
507687cd 770 TH2F* hPrimaryNSD = (TH2F*)listOfPrim->FindObject( Form("hPrimaryNSD_FMD_%c_vtx%d",ringChar,v));
771
ab3e0abc 772 TH2F* hAnalysed = (TH2F*)listOfPrim->FindObject(Form("Analysed_FMD%c_vtx%d",ringChar,v));
059c7c6b 773 TH2F* hAnalysedNSD = (TH2F*)listOfPrim->FindObject(Form("AnalysedNSD_FMD%c_vtx%d",ringChar,v));
ab3e0abc 774 TH2F* hInel = (TH2F*)listOfPrim->FindObject(Form("Inel_FMD%c_vtx%d",ringChar,v));
059c7c6b 775 TH2F* hNSD = (TH2F*)listOfPrim->FindObject(Form("NSD_FMD%c_vtx%d",ringChar,v));
d31dce33 776
ab3e0abc 777 fListOfPrimaries.Add(hPrimary);
df2a9c32 778 if(hPrimaryNSD)
779 fListOfPrimaries.Add(hPrimaryNSD);
ab3e0abc 780 fListOfPrimaries.Add(hAnalysed);
781 fListOfPrimaries.Add(hInel);
059c7c6b 782 fListOfPrimaries.Add(hAnalysedNSD);
783 fListOfPrimaries.Add(hNSD);
d31dce33 784 }
785 }
786 GenerateCorrection();
787
788 TFile fout("backgroundFromFile.root","recreate");
789 fListOfHits.Write();
790 fListOfPrimaries.Write();
791 fListOfCorrection.Write();
792 fVertexBins.Write();
793
794 fout.Close();
0b0a4ae5 795 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
d31dce33 796 if(storeInOCDB) {
0b0a4ae5 797 TFile fbg(pars->GetPath(pars->GetBackgroundID()),"RECREATE");
d31dce33 798 fBackground->Write(AliFMDAnaParameters::GetBackgroundID());
799 fbg.Close();
0b0a4ae5 800 TFile feselect(pars->GetPath(pars->GetEventSelectionEffID()),"RECREATE");
d31dce33 801 fEventSelectionEff->Write(AliFMDAnaParameters::GetEventSelectionEffID());
802 feselect.Close();
803
804 }
805
806
807
808
809}
810//_____________________________________________________________________
811//
812// EOF
813//