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