96957ca65b2ef6fff78f1f0e20e1e8a5b0308b30
[u/mrichter/AliRoot.git] / FMD / analysis / AliFMDAnalysisTaskDndeta.cxx
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>
9 #include "TH1F.h"
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"
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)
30
31
32 AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta()
33 : fDebug(0),
34   fOutputList(0),
35   fInputList(0),
36   fArray(0),
37   fInputArray(0),
38   fVertexString(0x0),
39   fNevents(),
40   fNMCevents(),
41   fStandalone(kTRUE),
42   fMCevent(0),
43   fLastTrackByStrip(),
44   fPrimary(kTRUE)
45 {
46   // Default constructor
47   DefineInput (0, TList::Class());
48   DefineOutput(0, TList::Class());
49 }
50 //_____________________________________________________________________
51 AliFMDAnalysisTaskDndeta::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(),
60     fNMCevents(),
61     fStandalone(kTRUE),
62     fMCevent(0),
63     fLastTrackByStrip(),
64     fPrimary(kTRUE)
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 //_____________________________________________________________________
75 void 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;
88   TH1F* hHits = 0;
89   TH1F* hPrimVertexBin = 0;
90   
91   
92   TH2F* hBgTmp   = pars->GetBackgroundCorrection(1, 'I', 0);
93   TH1F* hPrimary = new TH1F("hMultvsEta","hMultvsEta",
94                             hBgTmp->GetNbinsX(),
95                             hBgTmp->GetXaxis()->GetXmin(),
96                             hBgTmp->GetXaxis()->GetXmax());
97   hPrimary->Sumw2();
98   fOutputList->Add(hPrimary);
99   Int_t nVtxbins = pars->GetNvtxBins();
100   
101   
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());
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             
131             hMult->Sumw2();
132             hHits->Sumw2();
133             fOutputList->Add(hMult);
134             fOutputList->Add(hHits);
135             vtxArray->AddAtAndExpand(hMult,i);
136             
137           }
138         } 
139     }
140   
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),
145                               hBgTmp->GetNbinsX(),
146                               hBgTmp->GetXaxis()->GetXmin(),
147                               hBgTmp->GetXaxis()->GetXmax());
148     hPrimVertexBin->Sumw2();
149     fOutputList->Add(hPrimVertexBin);
150     
151   }
152   
153   fNevents.SetBins(nVtxbins,0,nVtxbins);
154   fNevents.SetName("nEvents");
155   fNMCevents.SetBins(nVtxbins,0,nVtxbins);
156   fNMCevents.SetName("nMCEvents");
157   fOutputList->Add(&fNevents);
158   fOutputList->Add(&fNMCevents);
159   
160 }
161 //_____________________________________________________________________
162 void AliFMDAnalysisTaskDndeta::ConnectInputData(Option_t */*option*/)
163 {
164   if(fStandalone) {
165     fInputList   = (TList*)GetInputData(0);
166     fVertexString = (TObjString*)GetInputData(1);
167   }
168 }
169 //_____________________________________________________________________
170 void AliFMDAnalysisTaskDndeta::Exec(Option_t */*option*/)
171 {
172   Int_t vtxbin   = fVertexString->GetString().Atoi();
173   fNevents.Fill(vtxbin);
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));
186       
187       hMultTotal->Add(hMultInput);
188       
189       
190     }
191   }
192   
193   if(fMCevent && fPrimary)
194     ProcessPrimary();
195   
196   if(fStandalone) {
197     PostData(0, fOutputList); 
198   }
199   
200 }
201 //_____________________________________________________________________
202 void AliFMDAnalysisTaskDndeta::Terminate(Option_t */*option*/) {
203   
204   
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);
214       Char_t ringChar = (ir == 0 ? 'I' : 'O');
215       for(Int_t i =0; i<nVtxbins; i++) {
216         TH2F* hMultTotal = (TH2F*)vtxArray->At(i);
217         TH1D* hMultProj   = hMultTotal->ProjectionX(Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ringChar,i),1,hMultTotal->GetNbinsY());
218         fOutputList->Add(hMultProj);
219       }
220     }
221   }
222   
223 }
224 //_____________________________________________________________________
225 void AliFMDAnalysisTaskDndeta::ProcessPrimary() {
226   
227   fLastTrackByStrip.Reset(-1);
228   
229   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
230   
231   AliMCParticle* particle = 0;
232   AliStack* stack = fMCevent->Stack();
233   
234   TH1F* hPrimary = (TH1F*)fOutputList->FindObject("hMultvsEta");
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;
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;
252     
253     if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
254       hPrimary->Fill(particle->Eta());
255       
256
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     }
264     
265     for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) {
266       
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 ) {
275         //Double_t x,y,z;
276         
277         Float_t   eta   = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));//-1*TMath::Log(TMath::Tan(0.5*theta));
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         
282         //if(det == 1 && ring == 'I')
283         //      std::cout<<"hit in "<<det<<"   "<<ring<<"   "<<sec<<"   "<<strip<<"   "<<std::endl;
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       }
293       }
294     
295     
296       }
297   
298 }
299 //_____________________________________________________________________
300 //
301 //
302 // EOF