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