Update to comply with the changes in AliFMDFloatMap
[u/mrichter/AliRoot.git] / FMD / analysis / AliFMDAnalysisTaskGenerateBackground.cxx
CommitLineData
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"
24ClassImp(AliFMDAnalysisTaskGenerateBackground)
25
26//_____________________________________________________________________
27AliFMDAnalysisTaskGenerateBackground::AliFMDAnalysisTaskGenerateBackground():
41bad769 28AliAnalysisTaskSE(),
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//_____________________________________________________________________
43AliFMDAnalysisTaskGenerateBackground::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//_____________________________________________________________________
63void 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//_____________________________________________________________________
114void AliFMDAnalysisTaskGenerateBackground::Init()
115{
2a667047 116 fLastTrackByStrip.Reset(-1);
117
118
119}
120//_____________________________________________________________________
121void 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//_____________________________________________________________________
217void 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//_____________________________________________________________________
229void 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 261void 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//