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