Various improvements
[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 #include <TFile.h>
29
30 ClassImp(AliFMDAnalysisTaskGenerateCorrection)
31 #if 0 
32 ; // This is for emacs - do not delete 
33 #endif
34
35 //_____________________________________________________________________
36 AliFMDAnalysisTaskGenerateCorrection::AliFMDAnalysisTaskGenerateCorrection():
37 AliAnalysisTaskSE(),
38   fListOfHits(), 
39   fListOfPrimaries(),
40   fListOfCorrection(),
41   fVertexBins(),
42   fLastTrackByStrip(0),
43   fHitsByStrip(0),
44   fZvtxCut(10),
45   fNvtxBins(10),
46   fNbinsEta(200),
47   fBackground(0),
48   fEventSelectionEff(0)
49 {
50   // Default constructor
51 }
52 //_____________________________________________________________________
53 AliFMDAnalysisTaskGenerateCorrection::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),
64   fBackground(0),
65   fEventSelectionEff(0)
66 {
67  
68   DefineOutput(1, TList::Class());
69   DefineOutput(2, TList::Class());
70   DefineOutput(3, TH1F::Class());
71   DefineOutput(4, TList::Class());
72 }
73 //_____________________________________________________________________
74 void AliFMDAnalysisTaskGenerateCorrection::UserCreateOutputObjects()
75 {
76 // Create the output containers
77 //
78   
79   std::cout<<"Creating output objects"<<std::endl;
80   
81   // SelectCollisionCandidates(AliVEvent::kAny);
82   
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),
87                                     fNbinsEta, -4,6, 20, 0,2*TMath::Pi());
88     hSPDhits->Sumw2();
89     fListOfHits.Add(hSPDhits);
90     
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     
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       
108       TH2F* hPrimary       = new TH2F(Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
109                                       Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
110                                       fNbinsEta, -4,6, nSec, 0,2*TMath::Pi());
111       hPrimary->Sumw2();
112       fListOfPrimaries.Add(hPrimary);
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       
119       
120       
121     }
122   }
123   
124  
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),
133                                   Form("DoubleHits_FMD%d%c",det,ringChar),fNbinsEta, -4,6);
134       TH1F* allHits = new TH1F(Form("allHits_FMD%d%c",det,ringChar),
135                                Form("allHits_FMD%d%c",det,ringChar), fNbinsEta, -4,6);
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),
145                                fNbinsEta, -4,6, nSec, 0, 2*TMath::Pi());
146         hHits->Sumw2();
147         fListOfHits.Add(hHits);
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         
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);
161   TH1F* hEventsAllNSD    = new TH1F("EventsAllNSD","EventsAllNSD",fNvtxBins,0,fNvtxBins);
162   TH1F* hEventsSelectedVtx  = new TH1F("EventsSelectedVtx","EventsSelectedVtx",fNvtxBins,0,fNvtxBins);  
163   TH1F* hEventsSelectedTrigger  = new TH1F("EventsSelectedTrigger","EventsSelectedTrigger",fNvtxBins,0,fNvtxBins);
164   
165   TH1F* hEventsSelectedNSDVtx  = new TH1F("EventsSelectedNSDVtx","EventsSelectedNSDVtx",fNvtxBins,0,fNvtxBins);  
166   TH1F* hEventsSelectedNSD     = new TH1F("EventsSelectedNSD","EventsSelectedNSD",fNvtxBins,0,fNvtxBins);
167   
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);
175   
176   hEventsSelected->Sumw2();
177   hEventsAll->Sumw2();
178   hEventsAllNSD->Sumw2();
179   
180   fListOfHits.Add(hEventsSelected);
181   fListOfHits.Add(hEventsSelectedVtx);
182   fListOfHits.Add(hEventsSelectedTrigger);
183   fListOfHits.Add(hEventsSelectedNSDVtx);
184   fListOfHits.Add(hEventsSelectedNSD);
185   
186   
187   fListOfPrimaries.Add(hEventsAll);
188   fListOfPrimaries.Add(hEventsAllNSD);
189   
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),
197                                        fNbinsEta, -4,6, nSec, 0,2*TMath::Pi());
198       
199       hAnalysed->Sumw2();
200       fListOfPrimaries.Add(hAnalysed);
201       
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       
209       TH2F* hInel           = new TH2F(Form("Inel_FMD%c_vtx%d",ringChar,v),
210                                        Form("Inel_FMD%c_vtx%d",ringChar,v),
211                                        fNbinsEta, -4,6, nSec, 0,2*TMath::Pi());
212       
213       hInel->Sumw2();
214       fListOfPrimaries.Add(hInel);
215       
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       
223     }
224   }
225   
226   fVertexBins.SetName("VertexBins");
227   fVertexBins.GetXaxis()->Set(fNvtxBins,-1*fZvtxCut,fZvtxCut);
228   
229  
230   
231 }
232 //_____________________________________________________________________
233 void AliFMDAnalysisTaskGenerateCorrection::Init()
234 {
235   fLastTrackByStrip.Reset(-1);
236   
237   
238 }
239 //_____________________________________________________________________
240 void 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();
251   
252   pars->SetTriggerStatus(esdevent);
253   
254   Double_t esdvertex[3];
255   Bool_t vtxStatus =  pars->GetVertex(esdevent,esdvertex);
256   
257   /* Double_t deltaEsd           = 2*fZvtxCut/fNvtxBins;
258   Double_t vertexBinDoubleEsd = (esdvertex[2] + fZvtxCut) / deltaEsd;
259   Int_t    vertexBinEsd       = (Int_t)vertexBinDoubleEsd;*/
260   
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   
271   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
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;
276     //return; 
277   }
278   Bool_t nsd = kTRUE;
279   if(pythiaGenHeader) {
280     Int_t pythiaType = pythiaGenHeader->ProcessType();
281     
282     if(pythiaType==92 || pythiaType==93)
283       nsd = kFALSE;
284   }
285   if(dpmHeader) {
286     Int_t processType = dpmHeader->ProcessType();
287     if(processType == 5 || processType == 6)  
288       nsd = kFALSE;
289     
290   }
291   
292   
293   TArrayF vertex;
294   genHeader->PrimaryVertex(vertex);
295   
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   
304   
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
311   TH1F* hEventsSelected           = (TH1F*)fListOfHits.FindObject("EventsSelected");
312   TH1F* hEventsSelectedVtx        = (TH1F*)fListOfHits.FindObject("EventsSelectedVtx");
313   TH1F* hEventsSelectedTrigger    = (TH1F*)fListOfHits.FindObject("EventsSelectedTrigger");
314   TH1F* hEventsSelectedNSDVtx     = (TH1F*)fListOfHits.FindObject("EventsSelectedNSDVtx");
315   TH1F* hEventsSelectedNSD        = (TH1F*)fListOfHits.FindObject("EventsSelectedNSD");
316   TH1F* hEventsAll                = (TH1F*)fListOfPrimaries.FindObject("EventsAll");
317   TH1F* hEventsAllNSD             = (TH1F*)fListOfPrimaries.FindObject("EventsAllNSD");
318   
319   // TH1F* hTriggered      = (TH1F*)fListOfHits.FindObject("Triggered");
320   //  TH1F* hTriggeredAll   = (TH1F*)fListOfPrimaries.FindObject("TriggeredAll");
321   
322   Bool_t vtxFound = kTRUE;
323   if(!vtxStatus)
324     vtxFound = kFALSE;
325   
326   //if(TMath::Abs(vertex.At(2)) > fZvtxCut) {
327   //  vtxFound = kFALSE;
328     
329   //}
330   
331   if(TMath::Abs(vertex.At(2)) > fZvtxCut) {
332     return;
333   }
334   
335   
336   Bool_t isTriggered    = pars->IsEventTriggered(AliFMDAnaParameters::kMB1);
337   Bool_t isTriggeredNSD = pars->IsEventTriggered(AliFMDAnaParameters::kNSD);
338   if(vtxFound && isTriggered) hEventsSelected->Fill(vertexBin);
339   
340   if(vtxFound) hEventsSelectedVtx->Fill(vertexBin);
341   if(isTriggered) hEventsSelectedTrigger->Fill(vertexBin);
342
343   if(vtxFound && isTriggeredNSD) hEventsSelectedNSDVtx->Fill(vertexBin);
344   if(isTriggeredNSD) hEventsSelectedNSD->Fill(vertexBin);
345
346   
347   hEventsAll->Fill(vertexBin);
348   if(nsd) hEventsAllNSD->Fill(vertexBin);
349   
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
376
377   for(Int_t i = 0 ;i<nTracks;i++) {
378     particle = (AliMCParticle*) mcevent->GetTrack(i);
379     
380     if(!particle)
381       continue;
382         
383     if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
384       
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         }
396       }
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));
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));
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));
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));
405       
406       //if(vtxFound && isTriggered) {
407       if(isTriggeredNSD) {
408         hAnalysedNSDInner->Fill(particle->Eta(),particle->Phi());
409         hAnalysedNSDOuter->Fill(particle->Eta(),particle->Phi());
410       }
411       if(isTriggered && vtxFound) {
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       
418       if(nsd) {
419         hNSDInner->Fill(particle->Eta(),particle->Phi());
420         hNSDOuter->Fill(particle->Eta(),particle->Phi());
421       }
422     }
423     if(!vtxFound || !isTriggered) continue;
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;
431
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);
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         
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   }
465   
466   
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));
474     
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   
480   PostData(1, &fListOfHits);
481   PostData(2, &fListOfPrimaries);
482   PostData(3, &fVertexBins);
483 }
484 //_____________________________________________________________________
485 void 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 //_____________________________________________________________________
497 void AliFMDAnalysisTaskGenerateCorrection::GenerateCorrection() {
498   
499   fBackground         = new AliFMDAnaCalibBackgroundCorrection();
500   fEventSelectionEff  = new AliFMDAnaCalibEventSelectionEfficiency();
501   
502   //Event selection
503   TH1F* hEventsSelected           = (TH1F*)fListOfHits.FindObject("EventsSelected");
504   TH1F* hEventsSelectedVtx        = (TH1F*)fListOfHits.FindObject("EventsSelectedVtx");
505   TH1F* hEventsSelectedTrigger    = (TH1F*)fListOfHits.FindObject("EventsSelectedTrigger");
506   TH1F* hEventsSelectedNSDVtx     = (TH1F*)fListOfHits.FindObject("EventsSelectedNSDVtx");
507   TH1F* hEventsSelectedNSD        = (TH1F*)fListOfHits.FindObject("EventsSelectedNSD");
508     
509   TH1F* hEventsAll                = (TH1F*)fListOfPrimaries.FindObject("EventsAll");
510   TH1F* hEventsAllNSD             = (TH1F*)fListOfPrimaries.FindObject("EventsAllNSD");
511   TH1F* hEventsSelectedVtxDivByTr = (TH1F*)hEventsSelectedVtx->Clone("hEventsSelectedVtxDivByTr");
512   TH1F* hEventsSelectedNSDVtxDivByNSD = (TH1F*)hEventsSelectedNSDVtx->Clone("hEventsSelectedNSDVtxDivByNSD");
513   
514   fListOfHits.Add(hEventsSelectedVtxDivByTr);
515   fListOfHits.Add(hEventsSelectedNSDVtxDivByNSD);
516   hEventsSelectedNSDVtxDivByNSD->Divide(hEventsSelectedNSD);
517   hEventsSelectedVtxDivByTr->Divide(hEventsSelectedTrigger);
518   
519   for(Int_t v=0;v<fNvtxBins ; v++) {
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));
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     
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));
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));
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))); 
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     
539     
540
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)));
543     
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     
547
548     hAnalysedInner->Divide(hInelInner);
549     hAnalysedOuter->Divide(hInelOuter);
550     
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     
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");
571     
572   }
573   
574   Float_t vtxEff = 1;
575   if(hEventsSelectedTrigger->GetEntries())
576     vtxEff = hEventsSelectedVtx->GetEntries() / hEventsSelectedTrigger->GetEntries();
577   fEventSelectionEff->SetVtxToTriggerRatio(vtxEff);
578   
579   //  hEventsAll->Divide(hEventsAll,hEventsSelected,1,1,"B");
580   hEventsSelectedVtx->Divide(hEventsAll);
581   hEventsSelectedTrigger->Divide(hEventsAll);
582   
583   hEventsSelectedNSDVtx->Divide(hEventsAllNSD);
584   hEventsSelectedNSD->Divide(hEventsAllNSD);
585   
586   for(Int_t i = 1; i<=hEventsSelected->GetNbinsX(); i++) {
587     if(hEventsSelected->GetBinContent(i) == 0 )
588       continue;
589     Float_t b    = hEventsSelected->GetBinContent(i);
590     // Float_t b    = hEventsSelected->GetBinContent(i)*hEventsSelectedTrigger->GetBinContent(i);
591     Float_t db   = hEventsSelected->GetBinError(i);
592     Float_t sum  = hEventsAll->GetBinContent(i);
593     Float_t dsum = hEventsAll->GetBinError(i);
594     Float_t a    = sum-b;
595     Float_t da   = TMath::Sqrt(TMath::Power(db,2) + TMath::Power(dsum,2));
596     
597     Float_t cor  = sum / b;
598     Float_t ecor = TMath::Sqrt(TMath::Power(da,2) + TMath::Power(a/(b*db),2)) / b;
599     
600     hEventsAll->SetBinContent(i,cor);
601     hEventsAll->SetBinError(i,ecor);
602     
603   }
604   
605   //TH1F* hEventTest = (TH1F*)hEventsAll->Clone("hEventTest");
606   
607   
608   
609   fEventSelectionEff->SetCorrection(hEventsAll);
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));
618       allHits->Divide(doubleHits);
619       
620       fBackground->SetDoubleHitCorrection(det,ring,allHits);
621       
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));
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));
626         hCorrection->Divide(hPrimary);
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));
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         }
635         
636         hCorrection->SetTitle(hCorrection->GetName());
637         fListOfCorrection.Add(hCorrection);
638         fBackground->SetBgCorrection(det,ring,vertexBin,hCorrection);
639         
640         if(hCorrectionNSD) {
641           fBackground->SetNSDBgCorrection(det,ring,vertexBin,hCorrectionNSD);
642           hCorrectionNSD->SetTitle(hCorrectionNSD->GetName());
643           fListOfCorrection.Add(hCorrectionNSD);
644         }
645         
646       }
647       
648     }
649   }
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));
653     if(!hSPDMult) continue;
654     
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);
659
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       
665     
666         
667       for(Int_t yy = 1; yy <=hSPDMult->GetNbinsY(); yy++) {
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         
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);
689   
690     fBackground->SetBgCorrection(0,'Q',vertexBin,hCorrection);
691   }
692   
693   
694   TAxis refAxis(fNvtxBins,-1*fZvtxCut,fZvtxCut);
695   fBackground->SetRefAxis(&refAxis);
696
697 }
698 //_____________________________________________________________________
699 void AliFMDAnalysisTaskGenerateCorrection::ReadFromFile(const Char_t* filename, Bool_t storeInOCDB, Int_t /*runNo*/) {
700   
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   
711   TH1F* hEventsSelected           = (TH1F*)listOfHits->FindObject("EventsSelected");
712   TH1F* hEventsSelectedVtx        = (TH1F*)listOfHits->FindObject("EventsSelectedVtx");
713   TH1F* hEventsSelectedTrigger    = (TH1F*)listOfHits->FindObject("EventsSelectedTrigger");
714   TH1F* hEventsSelectedNSDVtx     = (TH1F*)listOfHits->FindObject("EventsSelectedNSDVtx");
715   TH1F* hEventsSelectedNSD        = (TH1F*)listOfHits->FindObject("EventsSelectedNSD");  
716   TH1F* hEventsAll                = (TH1F*)listOfPrim->FindObject("EventsAll");
717   TH1F* hEventsAllNSD             = (TH1F*)listOfPrim->FindObject("EventsAllNSD");
718   
719   fListOfHits.Add(hEventsSelected);
720   fListOfHits.Add(hEventsSelectedVtx);
721   fListOfHits.Add(hEventsSelectedNSD);
722   fListOfHits.Add(hEventsSelectedNSDVtx);
723   fListOfHits.Add(hEventsSelectedTrigger);
724   fListOfPrimaries.Add(hEventsAll);
725   fListOfPrimaries.Add(hEventsAllNSD);
726   
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
734   for(Int_t det =1; det<=3;det++)
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);
749               TH2F* hHitsNSD       = (TH2F*)listOfHits->FindObject(Form("hHitsNSD_FMD%d%c_vtx%d", det,ringChar,v));
750               if(hHitsNSD)
751                 fListOfHits.Add(hHitsNSD);
752             }
753         }
754     }
755   for(Int_t v=0; v<fNvtxBins;v++) {
756     TH2F* hSPDHits          = (TH2F*)listOfHits->FindObject(Form("hSPDhits_vtx%d", v));   
757     fListOfHits.Add(hSPDHits);
758     TH1F* hFMDReadChannels  = (TH1F*)listOfHits->FindObject(Form("hFMDReadChannels_vtx%d",v));
759     TH1F* hFMDAllChannels   = (TH1F*)listOfHits->FindObject(Form("hFMDAllChannels_vtx%d",v));
760     
761     if(hFMDReadChannels && hFMDAllChannels) {
762       fListOfHits.Add(hFMDReadChannels);
763       fListOfHits.Add(hFMDAllChannels);
764     }
765     for(Int_t iring = 0; iring<2;iring++) {
766       Char_t ringChar = (iring == 0 ? 'I' : 'O');
767       
768       
769       TH2F* hPrimary       = (TH2F*)listOfPrim->FindObject( Form("hPrimary_FMD_%c_vtx%d",ringChar,v));
770       TH2F* hPrimaryNSD       = (TH2F*)listOfPrim->FindObject( Form("hPrimaryNSD_FMD_%c_vtx%d",ringChar,v));
771       
772       TH2F* hAnalysed      = (TH2F*)listOfPrim->FindObject(Form("Analysed_FMD%c_vtx%d",ringChar,v));
773       TH2F* hAnalysedNSD   = (TH2F*)listOfPrim->FindObject(Form("AnalysedNSD_FMD%c_vtx%d",ringChar,v));
774       TH2F* hInel          = (TH2F*)listOfPrim->FindObject(Form("Inel_FMD%c_vtx%d",ringChar,v));
775       TH2F* hNSD           = (TH2F*)listOfPrim->FindObject(Form("NSD_FMD%c_vtx%d",ringChar,v));
776       
777       fListOfPrimaries.Add(hPrimary);      
778       if(hPrimaryNSD)
779         fListOfPrimaries.Add(hPrimaryNSD); 
780       fListOfPrimaries.Add(hAnalysed);
781       fListOfPrimaries.Add(hInel);      
782       fListOfPrimaries.Add(hAnalysedNSD);
783       fListOfPrimaries.Add(hNSD);      
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();
795   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
796   if (!storeInOCDB) {
797     pars->SetBackgroundPath(".");
798     pars->SetEnergyPath(".");
799     pars->SetEventSelectionPath(".");
800     pars->SetSharingEfficiencyPath(".");
801   }
802   AliInfo(Form("Generating %s", pars->GetPath(pars->GetBackgroundID())));
803   TFile fbg(pars->GetPath(pars->GetBackgroundID()),"RECREATE");
804   fBackground->Write(AliFMDAnaParameters::GetBackgroundID());
805   fbg.Close();
806
807   AliInfo(Form("Generating %s", pars->GetPath(pars->GetEventSelectionEffID())));
808   TFile feselect(pars->GetPath(pars->GetEventSelectionEffID()),"RECREATE");
809   fEventSelectionEff->Write(AliFMDAnaParameters::GetEventSelectionEffID());
810   feselect.Close();
811 }
812 //_____________________________________________________________________
813 //
814 // EOF
815 //