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 "AliHeader.h"
26 //#include "TDatabasePDG.h"
27 //#include "TParticlePDG.h"
28 #include "AliFMDStripIndex.h"
29 ClassImp(AliFMDAnalysisTaskDndeta)
32 AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta()
42 // Default constructor
43 DefineInput (0, TList::Class());
44 DefineOutput(0, TList::Class());
46 //_____________________________________________________________________
47 AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta(const char* name, Bool_t SE):
48 AliAnalysisTask(name, "Density"),
60 DefineInput (0, TList::Class());
61 DefineInput(1, TObjString::Class());
62 DefineOutput(0, TList::Class());
66 //_____________________________________________________________________
67 void AliFMDAnalysisTaskDndeta::CreateOutputObjects()
69 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
72 fOutputList = new TList();
73 fOutputList->SetName("BackgroundCorrected");
78 TH1F* hPrimVertexBin = 0;
81 TH2F* hBgTmp = pars->GetBackgroundCorrection(1, 'I', 0);
82 TH1F* hPrimary = new TH1F("hMultvsEta","hMultvsEta",
84 hBgTmp->GetXaxis()->GetXmin(),
85 hBgTmp->GetXaxis()->GetXmax());
87 fOutputList->Add(hPrimary);
88 Int_t nVtxbins = pars->GetNvtxBins();
91 for(Int_t det =1; det<=3;det++)
93 Int_t nRings = (det==1 ? 1 : 2);
94 for(Int_t ring = 0;ring<nRings;ring++)
96 Char_t ringChar = (ring == 0 ? 'I' : 'O');
97 Int_t nSec = (ring == 0 ? 20 : 40);
100 for(Int_t i = 0; i< nVtxbins; i++) {
101 TH2F* hBg = pars->GetBackgroundCorrection(det, ringChar, i);
102 hMult = new TH2F(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),
104 hBg->GetXaxis()->GetXmin(),
105 hBg->GetXaxis()->GetXmax(),
106 nSec, 0, 2*TMath::Pi());
108 hHits = new TH1F(Form("hMCHits_FMD%d%c_vtxbin%d",det,ringChar,i),Form("hMCHits_FMD%d%c_vtxbin%d",det,ringChar,i),
110 hBg->GetXaxis()->GetXmin(),
111 hBg->GetXaxis()->GetXmax());
117 fOutputList->Add(hMult);
118 fOutputList->Add(hHits);
124 for(Int_t i = 0; i< nVtxbins; i++) {
126 hPrimVertexBin = new TH1F(Form("primmult_vtxbin%d",i),
127 Form("primmult_vtxbin%d",i),
129 hBgTmp->GetXaxis()->GetXmin(),
130 hBgTmp->GetXaxis()->GetXmax());
131 hPrimVertexBin->Sumw2();
132 fOutputList->Add(hPrimVertexBin);
136 fNevents.SetBins(nVtxbins,0,nVtxbins);
137 fNevents.SetName("nEvents");
138 fNMCevents.SetBins(nVtxbins,0,nVtxbins);
139 fNMCevents.SetName("nMCEvents");
140 fOutputList->Add(&fNevents);
141 fOutputList->Add(&fNMCevents);
144 //_____________________________________________________________________
145 void AliFMDAnalysisTaskDndeta::ConnectInputData(Option_t */*option*/)
148 fInputList = (TList*)GetInputData(0);
151 //_____________________________________________________________________
152 void AliFMDAnalysisTaskDndeta::Exec(Option_t */*option*/)
154 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
155 fVertexString = (TObjString*)fInputList->At(0);
157 Int_t vtxbin = fVertexString->GetString().Atoi();
159 fNevents.Fill(vtxbin);
161 for(UShort_t det=1;det<=3;det++) {
162 Int_t nRings = (det==1 ? 1 : 2);
163 for (UShort_t ir = 0; ir < nRings; ir++) {
164 Char_t ringChar = (ir == 0 ? 'I' : 'O');
166 TH2F* hMultTotal = (TH2F*)fOutputList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
169 TH2F* hMultInput = (TH2F*)fInputList->FindObject(Form("mult_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
171 hMultTotal->Add(hMultInput);
176 if(pars->GetProcessPrimary())
180 PostData(0, fOutputList);
184 //_____________________________________________________________________
185 void AliFMDAnalysisTaskDndeta::Terminate(Option_t */*option*/) {
188 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
190 Int_t nVtxbins = pars->GetNvtxBins();
192 for(UShort_t det=1;det<=3;det++) {
193 Int_t nRings = (det==1 ? 1 : 2);
194 for (UShort_t ir = 0; ir < nRings; ir++) {
195 Char_t ringChar = (ir == 0 ? 'I' : 'O');
196 for(Int_t i =0; i<nVtxbins; i++) {
198 TH2F* hMultTotal = (TH2F*)fOutputList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i));
199 TH2F* hMultTrVtx = (TH2F*)hMultTotal->Clone(Form("dNdeta_FMD%d%c_TrVtx_vtxbin%d",det,ringChar,i));
201 hMultTotal->Scale(pars->GetEventSelectionEfficiency(i));
202 TH1D* hMultProj = hMultTotal->ProjectionX(Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ringChar,i),1,hMultTotal->GetNbinsY());
203 TH1D* hMultProjTrVtx = hMultTrVtx->ProjectionX(Form("dNdeta_FMD%d%c_TrVtx_vtxbin%d_proj",det,ringChar,i),1,hMultTotal->GetNbinsY());
204 fOutputList->Add(hMultTrVtx);
205 fOutputList->Add(hMultProj);
206 fOutputList->Add(hMultProjTrVtx);
212 //_____________________________________________________________________
213 void AliFMDAnalysisTaskDndeta::ProcessPrimary() {
215 AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
216 AliMCEvent* mcEvent = eventHandler->MCEvent();
220 fLastTrackByStrip.Reset(-1);
222 AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
224 AliMCParticle* particle = 0;
225 AliStack* stack = mcEvent->Stack();
227 TH1F* hPrimary = (TH1F*)fOutputList->FindObject("hMultvsEta");
228 AliHeader* header = mcEvent->Header();
229 AliGenEventHeader* genHeader = header->GenEventHeader();
232 genHeader->PrimaryVertex(vertex);
233 if(TMath::Abs(vertex.At(2)) > pars->GetVtxCutZ())
235 Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
236 Double_t vertexBinDouble = (vertex.At(2) + pars->GetVtxCutZ()) / delta;
237 Int_t vertexBin = (Int_t)vertexBinDouble;
239 Bool_t firstTrack = kTRUE;
241 // we loop over the primaries only unless we need the hits (diagnostics running slowly)
242 Int_t nTracks = stack->GetNprimary();
243 if(pars->GetProcessHits())
244 nTracks = stack->GetNtrack();
246 for(Int_t i = 0 ;i<nTracks;i++) {
247 particle = mcEvent->GetTrack(i);
251 if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
252 hPrimary->Fill(particle->Eta());
255 TH1F* hPrimVtxBin = (TH1F*)fOutputList->FindObject(Form("primmult_vtxbin%d",vertexBin));
256 hPrimVtxBin->Fill(particle->Eta());
258 fNMCevents.Fill(vertexBin);
263 if(pars->GetProcessHits()) {
264 for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) {
266 AliTrackReference* ref = particle->GetTrackReference(j);
267 UShort_t det,sec,strip;
269 if(ref->DetectorId() != AliTrackReference::kFMD)
271 AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip);
272 Float_t thisStripTrack = fLastTrackByStrip(det,ring,sec,strip);
273 if(particle->Charge() != 0 && i != thisStripTrack ) {
276 Float_t eta = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));//-1*TMath::Log(TMath::Tan(0.5*theta));
277 TH1F* hHits = (TH1F*)fOutputList->FindObject(Form("hMCHits_FMD%d%c_vtxbin%d",det,ring,vertexBin));
279 Float_t nstrips = (ring =='O' ? 256 : 512);
281 fLastTrackByStrip(det,ring,sec,strip) = (Float_t)i;
284 fLastTrackByStrip(det,ring,sec,strip-1) = (Float_t)i;
285 if(strip < (nstrips - 1))
286 fLastTrackByStrip(det,ring,sec,strip+1) = (Float_t)i;
293 //_____________________________________________________________________