1 #include "AliFMDAnalysisTaskGenerateBackground.h"
2 #include "AliESDEvent.h"
6 #include "AliTrackReference.h"
8 #include "AliFMDAnaParameters.h"
9 #include "AliFMDStripIndex.h"
11 #include "AliMCParticle.h"
12 #include "AliMCEvent.h"
13 #include "AliFMDGeometry.h"
15 #include "AliGenEventHeader.h"
16 #include "AliHeader.h"
17 #include "AliFMDAnaCalibBackgroundCorrection.h"
18 #include "AliCDBManager.h"
20 #include "AliCDBMetaData.h"
24 ClassImp(AliFMDAnalysisTaskGenerateBackground)
26 //_____________________________________________________________________
27 AliFMDAnalysisTaskGenerateBackground::AliFMDAnalysisTaskGenerateBackground():
30 // Default constructor
32 //_____________________________________________________________________
33 AliFMDAnalysisTaskGenerateBackground::AliFMDAnalysisTaskGenerateBackground(const char* name):
34 AliAnalysisTaskSE(name),
40 DefineOutput(1, TList::Class());
41 DefineOutput(2, TList::Class());
42 DefineOutput(3, TH1F::Class());
43 DefineOutput(4, TList::Class());
45 //_____________________________________________________________________
46 void AliFMDAnalysisTaskGenerateBackground::UserCreateOutputObjects()
48 // Create the output containers
51 std::cout<<"Creating output objects"<<std::endl;
52 for(Int_t iring = 0; iring<2;iring++) {
53 Char_t ringChar = (iring == 0 ? 'I' : 'O');
54 Int_t nSec = (iring == 1 ? 40 : 20);
55 for(Int_t v=0; v<fNvtxBins;v++) {
57 TH2F* hPrimary = new TH2F(Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
58 Form("hPrimary_FMD_%c_vtx%d",ringChar,v),
59 fNbinsEta, -6,6, nSec, 0,2*TMath::Pi());
61 fListOfPrimaries.Add(hPrimary);
66 for(Int_t det =1; det<=3;det++) {
67 Int_t nRings = (det==1 ? 1 : 2);
68 for(Int_t ring = 0;ring<nRings;ring++) {
69 Int_t nSec = (ring == 1 ? 40 : 20);
70 Char_t ringChar = (ring == 0 ? 'I' : 'O');
71 for(Int_t v=0; v<fNvtxBins;v++) {
72 TH2F* hHits = new TH2F(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
73 Form("hHits_FMD%d%c_vtx%d", det,ringChar,v),
74 fNbinsEta, -6,6, nSec, 0, 2*TMath::Pi());
76 fListOfHits.Add(hHits);
81 TH1F* doubleHits = new TH1F("DoubleHits", "DoubleHits",
83 TH1F* allHits = new TH1F("allHits", "allHits",
85 fVertexBins.SetName("VertexBins");
86 fVertexBins.GetXaxis()->Set(fNvtxBins,-1*fZvtxCut,fZvtxCut);
87 fListOfHits.Add(allHits);
88 fListOfHits.Add(doubleHits);
90 //_____________________________________________________________________
91 void AliFMDAnalysisTaskGenerateBackground::Init()
93 std::cout<<"Init"<<std::endl;
96 fLastTrackByStrip.Reset(-1);
100 //_____________________________________________________________________
101 void AliFMDAnalysisTaskGenerateBackground::UserExec(Option_t */*option*/)
104 fLastTrackByStrip.Reset(-1);
105 fHitsByStrip.Reset(0);
106 AliMCEvent* mcevent = MCEvent();
107 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
109 AliMCParticle* particle = 0;
110 AliStack* stack = mcevent->Stack();
112 UShort_t det,sec,strip;
115 Int_t nTracks = mcevent->GetNumberOfTracks();
116 AliHeader* header = mcevent->Header();
117 AliGenEventHeader* genHeader = header->GenEventHeader();
120 genHeader->PrimaryVertex(vertex);
122 if(TMath::Abs(vertex.At(2)) > fZvtxCut)
124 for(Int_t i = 0 ;i<nTracks;i++) {
125 particle = mcevent->GetTrack(i);
130 Double_t delta = 2*fZvtxCut/fNvtxBins;
131 Double_t vertexBinDouble = (vertex.At(2) + fZvtxCut) / delta;
132 Int_t vertexBin = (Int_t)vertexBinDouble;
134 if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
137 TH2F* hPrimaryInner = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",'I',vertexBin));
138 TH2F* hPrimaryOuter = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",'O',vertexBin));
139 hPrimaryInner->Fill(particle->Eta(),particle->Phi());
140 hPrimaryOuter->Fill(particle->Eta(),particle->Phi());
143 for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) {
145 AliTrackReference* ref = particle->GetTrackReference(j);
147 if(ref->DetectorId() != AliTrackReference::kFMD)
149 AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip);
150 Float_t thisStripTrack = fLastTrackByStrip.operator()(det,ring,sec,strip);
151 if(particle->Charge() != 0 && i != thisStripTrack ) {
153 AliFMDGeometry* fmdgeo = AliFMDGeometry::Instance();
154 fmdgeo->Detector2XYZ(det,ring,sec,strip,x,y,z);
156 Float_t phi = TMath::ATan2(y,x);
157 if(phi<0) phi = phi+2*TMath::Pi();
158 Float_t r = TMath::Sqrt(TMath::Power(x,2)+TMath::Power(y,2));
159 Float_t theta = TMath::ATan2(r,z-vertex.At(2));
160 Float_t eta = -1*TMath::Log(TMath::Tan(0.5*theta));
161 TH2F* hHits = (TH2F*)fListOfHits.FindObject(Form("hHits_FMD%d%c_vtx%d", det,ring,vertexBin));
162 hHits->Fill(eta,phi);
163 Float_t nstrips = (ring =='O' ? 256 : 512);
164 fHitsByStrip.operator()(det,ring,sec,strip) +=1;
165 TH1F* allHits = (TH1F*)fListOfHits.FindObject("allHits");
166 TH1F* doubleHits = (TH1F*)fListOfHits.FindObject("DoubleHits");
168 if(fHitsByStrip.operator()(det,ring,sec,strip) == 1)
171 doubleHits->Fill(eta);
172 /*if(fHitsByStrip.operator()(det,ring,sec,strip) == 2){
173 TH1F* doubleHits = (TH1F*)fListOfHits.FindObject("DoubleHits");
174 doubleHits->Fill(eta,2);
176 //if(fHitsByStrip.operator()(det,ring,sec,strip) > 1){
177 // doubleHits->Fill(eta);
181 fLastTrackByStrip.operator()(det,ring,sec,strip) = (Float_t)i;
183 fLastTrackByStrip.operator()(det,ring,sec,strip-1) = (Float_t)i;
184 if(strip < (nstrips - 1))
185 fLastTrackByStrip.operator()(det,ring,sec,strip+1) = (Float_t)i;
192 PostData(1, &fListOfHits);
193 PostData(2, &fListOfPrimaries);
194 PostData(3, &fVertexBins);
196 //_____________________________________________________________________
197 void AliFMDAnalysisTaskGenerateBackground::Terminate(Option_t */*option*/)
199 TH1F* allHits = (TH1F*)fListOfHits.FindObject("allHits");
200 TH1F* doubleHits = (TH1F*)fListOfHits.FindObject("DoubleHits");
202 doubleHits->Divide(allHits);
203 GenerateCorrection();
204 PostData(1, &fListOfHits);
205 PostData(4, &fListOfCorrection);
208 //_____________________________________________________________________
209 void AliFMDAnalysisTaskGenerateBackground::GenerateCorrection() {
211 fBackground = new AliFMDAnaCalibBackgroundCorrection();
213 for(Int_t det= 1; det <=3; det++) {
214 Int_t nRings = (det==1 ? 1 : 2);
216 for(Int_t iring = 0; iring<nRings; iring++) {
217 Char_t ring = (iring == 0 ? 'I' : 'O');
218 for(Int_t vertexBin=0;vertexBin<fNvtxBins ;vertexBin++) {
219 TH2F* hHits = (TH2F*)fListOfHits.FindObject(Form("hHits_FMD%d%c_vtx%d", det,ring,vertexBin));
220 TH2F* hPrimary = (TH2F*)fListOfPrimaries.FindObject( Form("hPrimary_FMD_%c_vtx%d",ring,vertexBin));
221 TH2F* hCorrection = (TH2F*)hHits->Clone(Form("FMD%d%c_vtxbin_%d_correction",det,ring,vertexBin));
222 hCorrection->Divide(hPrimary);
223 hCorrection->SetTitle(hCorrection->GetName());
224 fListOfCorrection.Add(hCorrection);
225 fBackground->SetBgCorrection(det,ring,vertexBin,hCorrection);
230 TAxis refAxis(fNvtxBins,-1*fZvtxCut,fZvtxCut);
231 fBackground->SetRefAxis(&refAxis);
234 //_____________________________________________________________________
235 void AliFMDAnalysisTaskGenerateBackground::ReadFromFile(const Char_t* filename, Bool_t storeInOCDB, Int_t runNo) {
237 TFile infile(filename);
238 TH1F* hVertex = (TH1F*)infile.Get("VertexBins");
239 fZvtxCut = hVertex->GetXaxis()->GetXmax();
240 fNvtxBins = hVertex->GetXaxis()->GetNbins();
241 fVertexBins.SetName("VertexBins");
242 fVertexBins.GetXaxis()->Set(fNvtxBins,-1*fZvtxCut,fZvtxCut);
244 TList* listOfHits = (TList*)infile.Get("Hits");
245 TList* listOfPrim = (TList*)infile.Get("Primaries");
247 for(Int_t det =1; det<=3;det++)
249 Int_t nRings = (det==1 ? 1 : 2);
250 for(Int_t ring = 0;ring<nRings;ring++)
252 Char_t ringChar = (ring == 0 ? 'I' : 'O');
253 for(Int_t v=0; v<fNvtxBins;v++)
256 TH2F* hHits = (TH2F*)listOfHits->FindObject(Form("hHits_FMD%d%c_vtx%d", det,ringChar,v));
257 fListOfHits.Add(hHits);
261 for(Int_t iring = 0; iring<2;iring++) {
262 Char_t ringChar = (iring == 0 ? 'I' : 'O');
263 for(Int_t v=0; v<fNvtxBins;v++) {
265 TH2F* hPrimary = (TH2F*)listOfPrim->FindObject( Form("hPrimary_FMD_%c_vtx%d",ringChar,v));
266 fListOfPrimaries.Add(hPrimary);
270 GenerateCorrection();
272 TFile fout("backgroundFromFile.root","recreate");
274 fListOfPrimaries.Write();
275 fListOfCorrection.Write();
280 AliCDBManager* cdb = AliCDBManager::Instance();
281 cdb->SetDefaultStorage("local://$ALICE_ROOT/OCDB");
282 AliCDBId id(AliFMDAnaParameters::GetBackgroundPath(),runNo,999999999);
284 AliCDBMetaData* meta = new AliCDBMetaData;
285 meta->SetResponsible(gSystem->GetUserInfo()->fRealName.Data());
286 meta->SetAliRootVersion(gROOT->GetVersion());
287 meta->SetBeamPeriod(1);
288 meta->SetComment("Background Correction for FMD");
289 meta->SetProperty("key1", fBackground );
290 cdb->Put(fBackground, id, meta);
295 //_____________________________________________________________________