2a667047 |
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" |
78f6f750 |
13 | //#include "AliFMDGeometry.h" |
2a667047 |
14 | #include "TArray.h" |
15 | #include "AliGenEventHeader.h" |
16 | #include "AliHeader.h" |
17 | #include "AliFMDAnaCalibBackgroundCorrection.h" |
78f6f750 |
18 | //#include "AliCDBManager.h" |
19 | //#include "AliCDBId.h" |
20 | //#include "AliCDBMetaData.h" |
2a667047 |
21 | #include "TSystem.h" |
22 | #include "TROOT.h" |
23 | #include "TAxis.h" |
24 | ClassImp(AliFMDAnalysisTaskGenerateBackground) |
25 | |
26 | //_____________________________________________________________________ |
27 | AliFMDAnalysisTaskGenerateBackground::AliFMDAnalysisTaskGenerateBackground(): |
41bad769 |
28 | AliAnalysisTaskSE(), |
29 | fListOfHits(), |
30 | fListOfPrimaries(), |
31 | fListOfCorrection(), |
32 | fVertexBins(), |
b83cd7a7 |
33 | fLastTrackByStrip(0), |
34 | fHitsByStrip(0), |
41bad769 |
35 | fZvtxCut(10), |
36 | fNvtxBins(10), |
37 | fNbinsEta(200), |
38 | fBackground(0) |
2a667047 |
39 | { |
40 | // Default constructor |
41 | } |
42 | //_____________________________________________________________________ |
43 | AliFMDAnalysisTaskGenerateBackground::AliFMDAnalysisTaskGenerateBackground(const char* name): |
44 | AliAnalysisTaskSE(name), |
41bad769 |
45 | fListOfHits(), |
46 | fListOfPrimaries(), |
47 | fListOfCorrection(), |
48 | fVertexBins(), |
b83cd7a7 |
49 | fLastTrackByStrip(0), |
50 | fHitsByStrip(0), |
2a667047 |
51 | fZvtxCut(10), |
52 | fNvtxBins(10), |
41bad769 |
53 | fNbinsEta(200), |
54 | fBackground(0) |
2a667047 |
55 | { |
f58a4769 |
56 | |
2a667047 |
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'); |
4fb49e43 |
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 | |
2a667047 |
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 | } |
4fb49e43 |
108 | |
2a667047 |
109 | fVertexBins.SetName("VertexBins"); |
110 | fVertexBins.GetXaxis()->Set(fNvtxBins,-1*fZvtxCut,fZvtxCut); |
4fb49e43 |
111 | |
2a667047 |
112 | } |
113 | //_____________________________________________________________________ |
114 | void AliFMDAnalysisTaskGenerateBackground::Init() |
115 | { |
2a667047 |
116 | fLastTrackByStrip.Reset(-1); |
117 | |
118 | |
119 | } |
120 | //_____________________________________________________________________ |
121 | void AliFMDAnalysisTaskGenerateBackground::UserExec(Option_t */*option*/) |
122 | { |
123 | |
124 | fLastTrackByStrip.Reset(-1); |
f58a4769 |
125 | fHitsByStrip.Reset(0); |
2a667047 |
126 | AliMCEvent* mcevent = MCEvent(); |
f58a4769 |
127 | AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance(); |
2a667047 |
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); |
b83cd7a7 |
170 | Float_t thisStripTrack = fLastTrackByStrip(det,ring,sec,strip); |
2a667047 |
171 | if(particle->Charge() != 0 && i != thisStripTrack ) { |
78f6f750 |
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)); |
2a667047 |
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); |
b83cd7a7 |
185 | fHitsByStrip(det,ring,sec,strip) +=1; |
4fb49e43 |
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)); |
f58a4769 |
188 | |
b83cd7a7 |
189 | if(fHitsByStrip(det,ring,sec,strip) == 1) |
f58a4769 |
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 | |
b83cd7a7 |
201 | fLastTrackByStrip(det,ring,sec,strip) = (Float_t)i; |
2a667047 |
202 | if(strip >0) |
b83cd7a7 |
203 | fLastTrackByStrip(det,ring,sec,strip-1) = (Float_t)i; |
2a667047 |
204 | if(strip < (nstrips - 1)) |
b83cd7a7 |
205 | fLastTrackByStrip(det,ring,sec,strip+1) = (Float_t)i; |
2a667047 |
206 | } |
207 | } |
208 | |
209 | } |
210 | |
f58a4769 |
211 | |
2a667047 |
212 | PostData(1, &fListOfHits); |
213 | PostData(2, &fListOfPrimaries); |
214 | PostData(3, &fVertexBins); |
215 | } |
216 | //_____________________________________________________________________ |
217 | void AliFMDAnalysisTaskGenerateBackground::Terminate(Option_t */*option*/) |
218 | { |
78f6f750 |
219 | /* TH1F* allHits = (TH1F*)fListOfHits.FindObject("allHits"); |
f58a4769 |
220 | TH1F* doubleHits = (TH1F*)fListOfHits.FindObject("DoubleHits"); |
221 | |
222 | doubleHits->Divide(allHits); |
2a667047 |
223 | GenerateCorrection(); |
f58a4769 |
224 | PostData(1, &fListOfHits); |
78f6f750 |
225 | PostData(4, &fListOfCorrection);*/ |
f58a4769 |
226 | |
2a667047 |
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'); |
4fb49e43 |
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); |
2a667047 |
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); |
4fb49e43 |
250 | |
251 | |
2a667047 |
252 | } |
253 | |
254 | } |
255 | } |
256 | TAxis refAxis(fNvtxBins,-1*fZvtxCut,fZvtxCut); |
257 | fBackground->SetRefAxis(&refAxis); |
258 | |
259 | } |
260 | //_____________________________________________________________________ |
41bad769 |
261 | void AliFMDAnalysisTaskGenerateBackground::ReadFromFile(const Char_t* filename, Bool_t storeInOCDB, Int_t /*runNo*/) { |
2a667047 |
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'); |
4fb49e43 |
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); |
2a667047 |
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(); |
4fb49e43 |
307 | |
2a667047 |
308 | fout.Close(); |
309 | |
310 | if(storeInOCDB) { |
78f6f750 |
311 | TFile fcalib("$ALICE_ROOT/FMD/Correction/Background/background.root","RECREATE"); |
d05586f1 |
312 | fBackground->Write(AliFMDAnaParameters::GetBackgroundID()); |
78f6f750 |
313 | fcalib.Close(); |
314 | /* AliCDBManager* cdb = AliCDBManager::Instance(); |
2a667047 |
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); |
78f6f750 |
325 | */ |
2a667047 |
326 | } |
327 | |
328 | } |
329 | //_____________________________________________________________________ |
330 | // |
331 | // EOF |
332 | // |