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