4 #include <TInterpreter.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"
20 #include "AliESDVertex.h"
22 #include "AliFMDAnaParameters.h"
23 //#include "AliFMDGeometry.h"
24 #include "AliGenEventHeader.h"
25 #include "AliGenPythiaEventHeader.h"
26 #include "AliHeader.h"
27 //#include "TDatabasePDG.h"
28 //#include "TParticlePDG.h"
29 #include "AliFMDStripIndex.h"
30 #include "AliESDInputHandler.h"
31 #include "AliGenDPMjetEventHeader.h"
33 ClassImp(AliFMDAnalysisTaskDndeta)
36 AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta()
50 // Default constructor
51 DefineInput (0, TList::Class());
52 DefineOutput(0, TList::Class());
54 //_____________________________________________________________________
55 AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta(const char* name, Bool_t SE):
56 AliAnalysisTask(name, "Density"),
72 DefineInput (0, TList::Class());
73 DefineInput(1, TObjString::Class());
74 DefineOutput(0, TList::Class());
78 //_____________________________________________________________________
79 void AliFMDAnalysisTaskDndeta::CreateOutputObjects()
81 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
84 fOutputList = new TList();
85 fOutputList->SetName("BackgroundCorrected");
92 TH1F* hPrimVertexBin = 0;
93 TH1F* hPrimVertexBinNSD = 0;
95 TH2F* hBgTmp = pars->GetBackgroundCorrection(1, 'I', 0);
96 TH1F* hPrimary = new TH1F("hMultvsEta","hMultvsEta",
98 hBgTmp->GetXaxis()->GetXmin(),
99 hBgTmp->GetXaxis()->GetXmax());
101 fOutputList->Add(hPrimary);
102 TH1F* hPrimaryNSD = new TH1F("hMultvsEtaNSD","hMultvsEtaNSD",
104 hBgTmp->GetXaxis()->GetXmin(),
105 hBgTmp->GetXaxis()->GetXmax());
106 hPrimaryNSD->Sumw2();
107 fOutputList->Add(hPrimaryNSD);
108 Int_t nVtxbins = pars->GetNvtxBins();
110 for(Int_t i = 0; i< nVtxbins; i++) {
112 for(Int_t det =1; det<=3;det++)
114 Int_t nRings = (det==1 ? 1 : 2);
115 for(Int_t ring = 0;ring<nRings;ring++)
117 Char_t ringChar = (ring == 0 ? 'I' : 'O');
118 Int_t nSec = (ring == 0 ? 20 : 40);
122 hBg = pars->GetBackgroundCorrection(det, ringChar, i);
123 hMult = new TH2F(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),
125 hBg->GetXaxis()->GetXmin(),
126 hBg->GetXaxis()->GetXmax(),
127 nSec, 0, 2*TMath::Pi());
128 hMultTrVtx = new TH2F(Form("dNdetaTrVtx_FMD%d%c_vtxbin%d",det,ringChar,i),Form("dNdetaTrVtx_FMD%d%c_vtxbin%d",det,ringChar,i),
130 hBg->GetXaxis()->GetXmin(),
131 hBg->GetXaxis()->GetXmax(),
132 nSec, 0, 2*TMath::Pi());
133 hMultNSD = new TH2F(Form("dNdetaNSD_FMD%d%c_vtxbin%d",det,ringChar,i),Form("dNdetaNSD_FMD%d%c_vtxbin%d",det,ringChar,i),
135 hBg->GetXaxis()->GetXmin(),
136 hBg->GetXaxis()->GetXmax(),
137 nSec, 0, 2*TMath::Pi());
138 hHits = new TH1F(Form("hMCHits_FMD%d%c_vtxbin%d",det,ringChar,i),Form("hMCHits_FMD%d%c_vtxbin%d",det,ringChar,i),
140 hBg->GetXaxis()->GetXmin(),
141 hBg->GetXaxis()->GetXmax());
143 fOutputList->Add(hHits);
146 fOutputList->Add(hMult);
149 fOutputList->Add(hMultTrVtx);
152 fOutputList->Add(hMultNSD);
158 for(Int_t i = 0; i< nVtxbins; i++) {
160 hPrimVertexBin = new TH1F(Form("primmult_vtxbin%d",i),
161 Form("primmult_vtxbin%d",i),
163 hBgTmp->GetXaxis()->GetXmin(),
164 hBgTmp->GetXaxis()->GetXmax());
165 hPrimVertexBin->Sumw2();
166 fOutputList->Add(hPrimVertexBin);
168 hPrimVertexBinNSD = new TH1F(Form("primmult_NSD_vtxbin%d",i),
169 Form("primmult_NSD_vtxbin%d",i),
171 hBgTmp->GetXaxis()->GetXmin(),
172 hBgTmp->GetXaxis()->GetXmax());
173 hPrimVertexBinNSD->Sumw2();
174 fOutputList->Add(hPrimVertexBinNSD);
178 TH2F* hSPDMult = new TH2F(Form("dNdeta_SPD_vtxbin%d",i),Form("dNdeta_SPD_vtxbin%d",i),
180 hBgTmp->GetXaxis()->GetXmin(),
181 hBgTmp->GetXaxis()->GetXmax(),
182 20, 0, 2*TMath::Pi());
184 fOutputList->Add(hSPDMult);
185 TH2F* hSPDMultTrVtx = new TH2F(Form("dNdetaTrVtx_SPD_vtxbin%d",i),Form("dNdetaTrVtx_SPD_vtxbin%d",i),
187 hBgTmp->GetXaxis()->GetXmin(),
188 hBgTmp->GetXaxis()->GetXmax(),
189 20, 0, 2*TMath::Pi());
190 hSPDMultTrVtx->Sumw2();
191 fOutputList->Add(hSPDMultTrVtx);
193 TH2F* hSPDMultNSD = new TH2F(Form("dNdetaNSD_SPD_vtxbin%d",i),Form("dNdetaNSD_SPD_vtxbin%d",i),
195 hBgTmp->GetXaxis()->GetXmin(),
196 hBgTmp->GetXaxis()->GetXmax(),
197 20, 0, 2*TMath::Pi());
198 hSPDMultNSD->Sumw2();
199 fOutputList->Add(hSPDMultNSD);
203 //AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
205 // TH2F* dNdetadphiHistogramTotal = new TH2F("dNdetadphiHistogramTotal","dNdetadphiHistogram;#eta;#Phi",pars->GetNetaBins(),-6,6,20,0,2*TMath::Pi());
206 //dNdetadphiHistogramTotal->SetErrorOption("g");
207 //fOutputList->Add(dNdetadphiHistogramTotal);
211 fNevents.SetBins(nVtxbins,0,nVtxbins);
212 fNevents.SetName("nEvents");
213 fNNSDevents.SetBins(nVtxbins,0,nVtxbins);
214 fNNSDevents.SetName("nNSDEvents");
216 fNMCevents.SetBins(nVtxbins,0,nVtxbins);
217 fNMCevents.SetName("nMCEvents");
218 fNMCNSDevents.SetBins(nVtxbins,0,nVtxbins);
219 fNMCNSDevents.SetName("nMCNSDEvents");
220 fOutputList->Add(&fNevents);
221 fOutputList->Add(&fNNSDevents);
222 fOutputList->Add(&fNMCevents);
223 fOutputList->Add(&fNMCNSDevents);
225 //_____________________________________________________________________
226 void AliFMDAnalysisTaskDndeta::ConnectInputData(Option_t */*option*/)
229 fInputList = (TList*)GetInputData(0);
232 //_____________________________________________________________________
233 void AliFMDAnalysisTaskDndeta::Exec(Option_t */*option*/)
235 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
236 fVertexString = (TObjString*)fInputList->At(0);
238 Int_t vtxbin = fVertexString->GetString().Atoi();
240 fNevents.Fill(vtxbin);
242 //AliESDInputHandler* eventHandler = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
243 //AliESDEvent* esd = eventHandler->GetEvent();
244 Bool_t nsd = pars->IsEventTriggered(AliFMDAnaParameters::kNSD);
245 if(nsd) fNNSDevents.Fill(vtxbin);
247 for(UShort_t det=1;det<=3;det++) {
248 Int_t nRings = (det==1 ? 1 : 2);
249 for (UShort_t ir = 0; ir < nRings; ir++) {
250 Char_t ringChar = (ir == 0 ? 'I' : 'O');
252 TH2F* hMultTotal = (TH2F*)fOutputList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
253 TH2F* hMultTotalTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
254 TH2F* hMultTotalNSD = (TH2F*)fOutputList->FindObject(Form("dNdetaNSD_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
257 TH2F* hMultInput = (TH2F*)fInputList->FindObject(Form("mult_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
258 TH2F* hMultInputTrVtx = (TH2F*)fInputList->FindObject(Form("multTrVtx_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
259 TH2F* hMultInputNSD = (TH2F*)fInputList->FindObject(Form("multNSD_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
261 hMultTotal->Add(hMultInput);
262 hMultTotalTrVtx->Add(hMultInputTrVtx);
264 hMultTotalNSD->Add(hMultInputNSD);
270 TH2F* hMultSPDTotal = (TH2F*)fOutputList->FindObject(Form("dNdeta_SPD_vtxbin%d",vtxbin));
271 TH2F* hMultSPDTotalTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_SPD_vtxbin%d",vtxbin));
272 TH2F* hMultSPDTotalNSD = (TH2F*)fOutputList->FindObject(Form("dNdetaNSD_SPD_vtxbin%d",vtxbin));
273 TH2F* hMultSPDInput = (TH2F*)fInputList->FindObject(Form("mult_SPD_vtxbin%d",vtxbin));
274 TH2F* hMultSPDInputTrVtx = (TH2F*)fInputList->FindObject(Form("multTrVtx_SPD_vtxbin%d",vtxbin));
275 TH2F* hMultSPDInputNSD = (TH2F*)fInputList->FindObject(Form("multNSD_SPD_vtxbin%d",vtxbin));
277 hMultSPDTotal->Add(hMultSPDInput);
278 hMultSPDTotalTrVtx->Add(hMultSPDInputTrVtx);
280 hMultSPDTotalNSD->Add(hMultSPDInputNSD);
282 if(pars->GetProcessPrimary())
285 //TH2F* dNdetadphiHistogram = (TH2F*)fOutputList->FindObject("dNdetadphiHistogramSPDTrVtx");
287 // TH2F* dNdetadphiHistogramTotal = (TH2F*)fOutputList->FindObject("dNdetadphiHistogramTotal");
290 // dNdetadphiHistogramTotal->Add(dNdetadphiHistogram);
296 PostData(0, fOutputList);
300 //_____________________________________________________________________
301 void AliFMDAnalysisTaskDndeta::Terminate(Option_t */*option*/) {
303 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
305 Int_t nVtxbins = pars->GetNvtxBins();
307 for(UShort_t det=1;det<=3;det++) {
308 Int_t nRings = (det==1 ? 1 : 2);
309 for (UShort_t ir = 0; ir < nRings; ir++) {
310 Char_t ringChar = (ir == 0 ? 'I' : 'O');
311 for(Int_t i =0; i<nVtxbins; i++) {
313 TH2F* hMultTotal = (TH2F*)fOutputList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i));
314 fOutputList->Add(hMultTotal->Clone(Form("%s_orig", hMultTotal->GetName())));
316 hMultTotal->Scale(fVtxEff);
319 TH2F* hMultTotalNSD = (TH2F*)fOutputList->FindObject(Form("dNdetaNSD_FMD%d%c_vtxbin%d",det,ringChar,i));
320 fOutputList->Add(hMultTotalNSD->Clone(Form("%s_orig", hMultTotalNSD->GetName())));
322 hMultTotalNSD->Scale(fVtxEffNSD);
325 //TH2F* hMultTrVtx = (TH2F*)hMultTotal->Clone(Form("dNdeta_FMD%d%c_TrVtx_vtxbin%d",det,ringChar,i));
326 TH2F* hMultTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_FMD%d%c_vtxbin%d",det,ringChar,i));
328 TH1D* hMultProj = hMultTotal->ProjectionX(Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ringChar,i),1,hMultTotal->GetNbinsY());
329 TH1D* hMultProjTrVtx = hMultTrVtx->ProjectionX(Form("dNdeta_FMD%d%c_TrVtx_vtxbin%d_proj",det,ringChar,i),1,hMultTrVtx->GetNbinsY());
330 TH1D* hMultProjNSD = hMultTotalNSD->ProjectionX(Form("dNdetaNSD_FMD%d%c_vtxbin%d_proj",det,ringChar,i),1,hMultTotalNSD->GetNbinsY());
332 //fOutputList->Add(hMultTrVtx);
333 fOutputList->Add(hMultProj);
334 fOutputList->Add(hMultProjTrVtx);
335 fOutputList->Add(hMultProjNSD);
340 for(Int_t i =0; i<nVtxbins; i++) {
342 TH2F* hSPDMult = (TH2F*)fOutputList->FindObject(Form("dNdeta_SPD_vtxbin%d",i));
343 TH2F* hSPDMultTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_SPD_vtxbin%d",i));
344 TH2F* hSPDMultNSD = (TH2F*)fOutputList->FindObject(Form("dNdetaNSD_SPD_vtxbin%d",i));
347 hSPDMult->Scale(fVtxEff);
350 hSPDMultNSD->Scale(fVtxEffNSD);
352 TH1D* hMultProj = hSPDMult->ProjectionX(Form("dNdeta_SPD_vtxbin%d_proj",i),1,hSPDMult->GetNbinsY());
353 TH1D* hMultProjTrVtx = hSPDMultTrVtx->ProjectionX(Form("dNdetaTrVtx_SPD_vtxbin%d_proj",i),1,hSPDMultTrVtx->GetNbinsY());
354 TH1D* hMultProjNSD = hSPDMultNSD->ProjectionX(Form("dNdetaNSD_SPD_vtxbin%d_proj",i),1,hSPDMultNSD->GetNbinsY());
355 fOutputList->Add(hMultProj);
356 fOutputList->Add(hMultProjTrVtx);
357 fOutputList->Add(hMultProjNSD);
361 std::cout<<"FMD analysis accepted "<<fNevents.GetEntries()<<" INEL events and "<<fNNSDevents.GetEntries()<<" NSD events"<<std::endl;
363 //_____________________________________________________________________
364 void AliFMDAnalysisTaskDndeta::ProcessPrimary() {
366 AliMCEventHandler* eventHandler =
367 dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()
368 ->GetMCtruthEventHandler());
369 if (!eventHandler) return;
371 AliMCEvent* mcEvent = eventHandler->MCEvent();
374 fLastTrackByStrip.Reset(-1);
376 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
378 AliMCParticle* particle = 0;
379 AliStack* stack = mcEvent->Stack();
381 TH1F* hPrimary = (TH1F*)fOutputList->FindObject("hMultvsEta");
382 TH1F* hPrimaryNSD = (TH1F*)fOutputList->FindObject("hMultvsEtaNSD");
383 AliHeader* header = mcEvent->Header();
384 AliGenEventHeader* genHeader = header->GenEventHeader();
386 AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
387 AliGenDPMjetEventHeader* dpmHeader = dynamic_cast<AliGenDPMjetEventHeader*>(header->GenEventHeader());
390 if (!pythiaGenHeader && !dpmHeader) {
391 std::cout<<" no pythia or dpm header!"<<std::endl;
394 if(pythiaGenHeader) {
395 Int_t pythiaType = pythiaGenHeader->ProcessType();
397 if(pythiaType==92||pythiaType==93)
402 Int_t processType = dpmHeader->ProcessType();
403 if(processType == 5 || processType == 6)
410 genHeader->PrimaryVertex(vertex);
411 if(TMath::Abs(vertex.At(2)) > pars->GetVtxCutZ())
413 Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
414 Double_t vertexBinDouble = (vertex.At(2) + pars->GetVtxCutZ()) / delta;
415 Int_t vertexBin = (Int_t)vertexBinDouble;
417 Bool_t firstTrack = kTRUE;
418 Bool_t firstTrackNSD = kTRUE;
420 TH1F* hPrimVtxBinNSD = (TH1F*)fOutputList->FindObject(Form("primmult_NSD_vtxbin%d",vertexBin));
421 TH1F* hPrimVtxBin = (TH1F*)fOutputList->FindObject(Form("primmult_vtxbin%d",vertexBin));
423 // we loop over the primaries only unless we need the hits (diagnostics running slowly)
424 Int_t nTracks = stack->GetNprimary();
425 if(pars->GetProcessHits())
426 nTracks = stack->GetNtrack();
428 for(Int_t i = 0 ;i<nTracks;i++) {
429 particle = (AliMCParticle*) mcEvent->GetTrack(i);
433 if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
434 hPrimary->Fill(particle->Eta());
435 hPrimVtxBin->Fill(particle->Eta());
437 fNMCevents.Fill(vertexBin);
442 hPrimaryNSD->Fill(particle->Eta());
443 hPrimVtxBinNSD->Fill(particle->Eta());
445 fNMCNSDevents.Fill(vertexBin);
446 firstTrackNSD = kFALSE;
453 if(pars->GetProcessHits()) {
455 for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) {
457 AliTrackReference* ref = particle->GetTrackReference(j);
458 UShort_t det,sec,strip;
460 if(ref->DetectorId() != AliTrackReference::kFMD)
462 AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip);
463 Float_t thisStripTrack = fLastTrackByStrip(det,ring,sec,strip);
464 if(particle->Charge() != 0 && i != thisStripTrack ) {
467 Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));//-1*TMath::Log(TMath::Tan(0.5*theta));
468 TH1F* hHits = (TH1F*)fOutputList->FindObject(Form("hMCHits_FMD%d%c_vtxbin%d",det,ring,vertexBin));
473 Float_t nstrips = (ring =='O' ? 256 : 512);
475 fLastTrackByStrip(det,ring,sec,strip) = (Float_t)i;
478 fLastTrackByStrip(det,ring,sec,strip-1) = (Float_t)i;
479 if(strip < (nstrips - 1))
480 fLastTrackByStrip(det,ring,sec,strip+1) = (Float_t)i;
494 //_____________________________________________________________________