bug fix: call of not yet implemented function removed
[u/mrichter/AliRoot.git] / FMD / analysis / AliFMDAnalysisTaskDndeta.cxx
CommitLineData
7c3e5162 1
2#include <TROOT.h>
3#include <TSystem.h>
4#include <TInterpreter.h>
5#include <TChain.h>
6#include <TFile.h>
7#include <TList.h>
8#include <iostream>
bb8a464f 9#include "TH1F.h"
7c3e5162 10#include "TH2F.h"
11#include "AliFMDAnalysisTaskDndeta.h"
12#include "AliAnalysisManager.h"
13#include "AliESDFMD.h"
14#include "AliESDEvent.h"
15#include "AliAODEvent.h"
16#include "AliAODHandler.h"
17#include "AliMCEventHandler.h"
18#include "AliStack.h"
19#include "AliLog.h"
20#include "AliESDVertex.h"
21#include "TMath.h"
22#include "AliFMDAnaParameters.h"
78f6f750 23//#include "AliFMDGeometry.h"
bb8a464f 24#include "AliGenEventHeader.h"
cc066cb9 25#include "AliHeader.h"
78f6f750 26//#include "TDatabasePDG.h"
27//#include "TParticlePDG.h"
cc066cb9 28#include "AliFMDStripIndex.h"
7c3e5162 29ClassImp(AliFMDAnalysisTaskDndeta)
30
31
32AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta()
33: fDebug(0),
34 fOutputList(0),
35 fInputList(0),
36 fArray(0),
37 fInputArray(0),
38 fVertexString(0x0),
39 fNevents(),
bb8a464f 40 fNMCevents(),
41 fStandalone(kTRUE),
cc066cb9 42 fMCevent(0),
b3546e91 43 fLastTrackByStrip(),
44 fPrimary(kTRUE)
7c3e5162 45{
46 // Default constructor
47 DefineInput (0, TList::Class());
48 DefineOutput(0, TList::Class());
49}
50//_____________________________________________________________________
51AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta(const char* name, Bool_t SE):
52 AliAnalysisTask(name, "Density"),
53 fDebug(0),
54 fOutputList(0),
55 fInputList(0),
56 fArray(),
57 fInputArray(0),
58 fVertexString(0x0),
59 fNevents(),
bb8a464f 60 fNMCevents(),
61 fStandalone(kTRUE),
cc066cb9 62 fMCevent(0),
b3546e91 63 fLastTrackByStrip(),
64 fPrimary(kTRUE)
7c3e5162 65{
66 fStandalone = SE;
67 if(fStandalone) {
68 DefineInput (0, TList::Class());
69 DefineInput(1, TObjString::Class());
70 DefineOutput(0, TList::Class());
71
72 }
73}
74//_____________________________________________________________________
75void AliFMDAnalysisTaskDndeta::CreateOutputObjects()
76{
77 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
78
79 fArray.SetName("FMD");
80 fArray.SetOwner();
81
82 if(!fOutputList)
83 fOutputList = new TList();
84 fOutputList->SetName("BackgroundCorrected");
85
86
87 TH2F* hMult = 0;
cc066cb9 88 TH1F* hHits = 0;
bb8a464f 89 TH1F* hPrimVertexBin = 0;
7c3e5162 90
bb8a464f 91
5ea1e0a9 92 TH2F* hBgTmp = pars->GetBackgroundCorrection(1, 'I', 0);
bb8a464f 93 TH1F* hPrimary = new TH1F("hMultvsEta","hMultvsEta",
5ea1e0a9 94 hBgTmp->GetNbinsX(),
95 hBgTmp->GetXaxis()->GetXmin(),
96 hBgTmp->GetXaxis()->GetXmax());
bb8a464f 97 hPrimary->Sumw2();
98 fOutputList->Add(hPrimary);
7c3e5162 99 Int_t nVtxbins = pars->GetNvtxBins();
100
bb8a464f 101
7c3e5162 102 for(Int_t det =1; det<=3;det++)
103 {
104 TObjArray* detArray = new TObjArray();
105 detArray->SetName(Form("FMD%d",det));
106 fArray.AddAtAndExpand(detArray,det);
107 Int_t nRings = (det==1 ? 1 : 2);
108 for(Int_t ring = 0;ring<nRings;ring++)
109 {
110 Char_t ringChar = (ring == 0 ? 'I' : 'O');
111 Int_t nSec = (ring == 0 ? 20 : 40);
112
113 TObjArray* vtxArray = new TObjArray();
114 vtxArray->SetName(Form("FMD%d%c",det,ringChar));
115 detArray->AddAtAndExpand(vtxArray,ring);
116 for(Int_t i = 0; i< nVtxbins; i++) {
117 TH2F* hBg = pars->GetBackgroundCorrection(det, ringChar, i);
118 hMult = new TH2F(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),
119 hBg->GetNbinsX(),
120 hBg->GetXaxis()->GetXmin(),
121 hBg->GetXaxis()->GetXmax(),
122 nSec, 0, 2*TMath::Pi());
cc066cb9 123
124 hHits = new TH1F(Form("hHits_FMD%d%c_vtxbin%d",det,ringChar,i),Form("hHits_FMD%d%c_vtxbin%d",det,ringChar,i),
125 hBg->GetNbinsX(),
126 hBg->GetXaxis()->GetXmin(),
127 hBg->GetXaxis()->GetXmax());
128
129
130
7c3e5162 131 hMult->Sumw2();
cc066cb9 132 hHits->Sumw2();
7c3e5162 133 fOutputList->Add(hMult);
cc066cb9 134 fOutputList->Add(hHits);
7c3e5162 135 vtxArray->AddAtAndExpand(hMult,i);
136
137 }
138 }
139 }
140
bb8a464f 141 for(Int_t i = 0; i< nVtxbins; i++) {
142
143 hPrimVertexBin = new TH1F(Form("primmult_vtxbin%d",i),
144 Form("primmult_vtxbin%d",i),
5ea1e0a9 145 hBgTmp->GetNbinsX(),
146 hBgTmp->GetXaxis()->GetXmin(),
147 hBgTmp->GetXaxis()->GetXmax());
bb8a464f 148 hPrimVertexBin->Sumw2();
149 fOutputList->Add(hPrimVertexBin);
150
151 }
152
7c3e5162 153 fNevents.SetBins(nVtxbins,0,nVtxbins);
154 fNevents.SetName("nEvents");
bb8a464f 155 fNMCevents.SetBins(nVtxbins,0,nVtxbins);
156 fNMCevents.SetName("nMCEvents");
7c3e5162 157 fOutputList->Add(&fNevents);
bb8a464f 158 fOutputList->Add(&fNMCevents);
7c3e5162 159
160}
161//_____________________________________________________________________
162void AliFMDAnalysisTaskDndeta::ConnectInputData(Option_t */*option*/)
163{
164 if(fStandalone) {
165 fInputList = (TList*)GetInputData(0);
166 fVertexString = (TObjString*)GetInputData(1);
167 }
168}
169//_____________________________________________________________________
170void AliFMDAnalysisTaskDndeta::Exec(Option_t */*option*/)
171{
172 Int_t vtxbin = fVertexString->GetString().Atoi();
173 fNevents.Fill(vtxbin);
7c3e5162 174 for(UShort_t det=1;det<=3;det++) {
175 //TObjArray* detInputArray = (TObjArray*)fInputArray->At(det);
176 TObjArray* detArray = (TObjArray*)fArray.At(det);
177 Int_t nRings = (det==1 ? 1 : 2);
178 for (UShort_t ir = 0; ir < nRings; ir++) {
179 Char_t ringChar = (ir == 0 ? 'I' : 'O');
180 //TObjArray* vtxInputArray = (TObjArray*)detInputArray->At(ir);
181 TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
182 TH2F* hMultTotal = (TH2F*)vtxArray->At(vtxbin);
183
184
185 TH2F* hMultInput = (TH2F*)fInputList->FindObject(Form("mult_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
1282ce49 186
7c3e5162 187 hMultTotal->Add(hMultInput);
188
189
190 }
191 }
bb8a464f 192
b3546e91 193 if(fMCevent && fPrimary)
bb8a464f 194 ProcessPrimary();
195
7c3e5162 196 if(fStandalone) {
197 PostData(0, fOutputList);
198 }
199
200}
201//_____________________________________________________________________
202void AliFMDAnalysisTaskDndeta::Terminate(Option_t */*option*/) {
203
1282ce49 204
7c3e5162 205 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
206
207 Int_t nVtxbins = pars->GetNvtxBins();
208
209 for(UShort_t det=1;det<=3;det++) {
210 TObjArray* detArray = (TObjArray*)fArray.At(det);
211 Int_t nRings = (det==1 ? 1 : 2);
212 for (UShort_t ir = 0; ir < nRings; ir++) {
213 TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
1282ce49 214 Char_t ringChar = (ir == 0 ? 'I' : 'O');
7c3e5162 215 for(Int_t i =0; i<nVtxbins; i++) {
216 TH2F* hMultTotal = (TH2F*)vtxArray->At(i);
1282ce49 217 TH1D* hMultProj = hMultTotal->ProjectionX(Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ringChar,i),1,hMultTotal->GetNbinsY());
218 fOutputList->Add(hMultProj);
7c3e5162 219 }
220 }
221 }
1282ce49 222
bb8a464f 223}
224//_____________________________________________________________________
225void AliFMDAnalysisTaskDndeta::ProcessPrimary() {
cc066cb9 226
227 fLastTrackByStrip.Reset(-1);
228
bb8a464f 229 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
230
231 AliMCParticle* particle = 0;
232 AliStack* stack = fMCevent->Stack();
233
234 TH1F* hPrimary = (TH1F*)fOutputList->FindObject("hMultvsEta");
cc066cb9 235 AliHeader* header = fMCevent->Header();
236 AliGenEventHeader* genHeader = header->GenEventHeader();
237
238 TArrayF vertex;
239 genHeader->PrimaryVertex(vertex);
240 if(TMath::Abs(vertex.At(2)) > pars->GetVtxCutZ())
241 return;
242 Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
243 Double_t vertexBinDouble = (vertex.At(2) + pars->GetVtxCutZ()) / delta;
244 Int_t vertexBin = (Int_t)vertexBinDouble;
bb8a464f 245
246 Bool_t firstTrack = kTRUE;
247 Int_t nTracks = fMCevent->GetNumberOfTracks();
248 for(Int_t i = 0 ;i<nTracks;i++) {
249 particle = fMCevent->GetTrack(i);
250 if(!particle)
251 continue;
6289b3e8 252
bb8a464f 253 if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
254 hPrimary->Fill(particle->Eta());
cc066cb9 255
6289b3e8 256
bb8a464f 257 TH1F* hPrimVtxBin = (TH1F*)fOutputList->FindObject(Form("primmult_vtxbin%d",vertexBin));
258 hPrimVtxBin->Fill(particle->Eta());
259 if(firstTrack) {
260 fNMCevents.Fill(vertexBin);
261 firstTrack = kFALSE;
262 }
263 }
b3546e91 264
cc066cb9 265 for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) {
bb8a464f 266
cc066cb9 267 AliTrackReference* ref = particle->GetTrackReference(j);
268 UShort_t det,sec,strip;
269 Char_t ring;
270 if(ref->DetectorId() != AliTrackReference::kFMD)
271 continue;
272 AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip);
273 Float_t thisStripTrack = fLastTrackByStrip.operator()(det,ring,sec,strip);
274 if(particle->Charge() != 0 && i != thisStripTrack ) {
78f6f750 275 //Double_t x,y,z;
cc066cb9 276
f58a4769 277 Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));//-1*TMath::Log(TMath::Tan(0.5*theta));
cc066cb9 278 TH1F* hHits = (TH1F*)fOutputList->FindObject(Form("hHits_FMD%d%c_vtxbin%d",det,ring,vertexBin));
279 hHits->Fill(eta);
280 Float_t nstrips = (ring =='O' ? 256 : 512);
281
f58a4769 282 //if(det == 1 && ring == 'I')
283 // std::cout<<"hit in "<<det<<" "<<ring<<" "<<sec<<" "<<strip<<" "<<std::endl;
cc066cb9 284 fLastTrackByStrip.operator()(det,ring,sec,strip) = (Float_t)i;
285
286 if(strip >0)
287 fLastTrackByStrip.operator()(det,ring,sec,strip-1) = (Float_t)i;
288 if(strip < (nstrips - 1))
289 fLastTrackByStrip.operator()(det,ring,sec,strip+1) = (Float_t)i;
290
291
292 }
b3546e91 293 }
cc066cb9 294
295
b3546e91 296 }
bb8a464f 297
7c3e5162 298}
299//_____________________________________________________________________
300//
301//
302// EOF