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