fixed the reading of the FMD primary data (speed issue)
[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(0),
44   fPrimary(kTRUE),
45   fRecordHits(kFALSE)
46 {
47   // Default constructor
48   DefineInput (0, TList::Class());
49   DefineOutput(0, TList::Class());
50 }
51 //_____________________________________________________________________
52 AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta(const char* name, Bool_t SE):
53     AliAnalysisTask(name, "Density"),
54     fDebug(0),
55     fOutputList(0),
56     fInputList(0),
57     fArray(),
58     fInputArray(0),
59     fVertexString(0x0),
60     fNevents(),
61     fNMCevents(),
62     fStandalone(kTRUE),
63     fMCevent(0),
64     fLastTrackByStrip(0),
65     fPrimary(kTRUE),
66     fRecordHits(kFALSE)
67 {
68   fStandalone = SE;
69   if(fStandalone) {
70     DefineInput (0, TList::Class());
71     DefineInput(1, TObjString::Class());
72     DefineOutput(0, TList::Class());
73     
74   }
75 }
76 //_____________________________________________________________________
77 void AliFMDAnalysisTaskDndeta::CreateOutputObjects()
78 {
79   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
80   
81   fArray.SetName("FMD");
82   fArray.SetOwner();
83   
84   if(!fOutputList)
85     fOutputList = new TList();
86   fOutputList->SetName("BackgroundCorrected");
87   
88   
89   TH2F* hMult = 0;
90   TH1F* hHits = 0;
91   TH1F* hPrimVertexBin = 0;
92   
93   
94   TH2F* hBgTmp   = pars->GetBackgroundCorrection(1, 'I', 0);
95   TH1F* hPrimary = new TH1F("hMultvsEta","hMultvsEta",
96                             hBgTmp->GetNbinsX(),
97                             hBgTmp->GetXaxis()->GetXmin(),
98                             hBgTmp->GetXaxis()->GetXmax());
99   hPrimary->Sumw2();
100   fOutputList->Add(hPrimary);
101   Int_t nVtxbins = pars->GetNvtxBins();
102   
103   
104   for(Int_t det =1; det<=3;det++)
105     {
106       TObjArray* detArray = new TObjArray();
107       detArray->SetName(Form("FMD%d",det));
108       fArray.AddAtAndExpand(detArray,det);
109       Int_t nRings = (det==1 ? 1 : 2);
110       for(Int_t ring = 0;ring<nRings;ring++)
111         {
112           Char_t ringChar = (ring == 0 ? 'I' : 'O');
113           Int_t  nSec     = (ring == 0 ? 20 : 40);
114           
115           TObjArray* vtxArray = new TObjArray();
116           vtxArray->SetName(Form("FMD%d%c",det,ringChar));
117           detArray->AddAtAndExpand(vtxArray,ring);
118           for(Int_t i = 0; i< nVtxbins; i++) {
119             TH2F* hBg = pars->GetBackgroundCorrection(det, ringChar, i);
120             hMult  = new TH2F(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),
121                               hBg->GetNbinsX(),
122                               hBg->GetXaxis()->GetXmin(),
123                               hBg->GetXaxis()->GetXmax(),
124                               nSec, 0, 2*TMath::Pi());
125             
126             hHits  = new TH1F(Form("hHits_FMD%d%c_vtxbin%d",det,ringChar,i),Form("hHits_FMD%d%c_vtxbin%d",det,ringChar,i),
127                               hBg->GetNbinsX(),
128                               hBg->GetXaxis()->GetXmin(),
129                               hBg->GetXaxis()->GetXmax());
130                             
131             
132             
133             hMult->Sumw2();
134             hHits->Sumw2();
135             fOutputList->Add(hMult);
136             fOutputList->Add(hHits);
137             vtxArray->AddAtAndExpand(hMult,i);
138             
139           }
140         } 
141     }
142   
143   for(Int_t i = 0; i< nVtxbins; i++) {
144    
145     hPrimVertexBin = new TH1F(Form("primmult_vtxbin%d",i),
146                               Form("primmult_vtxbin%d",i),
147                               hBgTmp->GetNbinsX(),
148                               hBgTmp->GetXaxis()->GetXmin(),
149                               hBgTmp->GetXaxis()->GetXmax());
150     hPrimVertexBin->Sumw2();
151     fOutputList->Add(hPrimVertexBin);
152     
153   }
154   
155   fNevents.SetBins(nVtxbins,0,nVtxbins);
156   fNevents.SetName("nEvents");
157   fNMCevents.SetBins(nVtxbins,0,nVtxbins);
158   fNMCevents.SetName("nMCEvents");
159   fOutputList->Add(&fNevents);
160   fOutputList->Add(&fNMCevents);
161   
162 }
163 //_____________________________________________________________________
164 void AliFMDAnalysisTaskDndeta::ConnectInputData(Option_t */*option*/)
165 {
166   if(fStandalone) {
167     fInputList   = (TList*)GetInputData(0);
168     fVertexString = (TObjString*)GetInputData(1);
169   }
170 }
171 //_____________________________________________________________________
172 void AliFMDAnalysisTaskDndeta::Exec(Option_t */*option*/)
173 {
174   Int_t vtxbin   = fVertexString->GetString().Atoi();
175   fNevents.Fill(vtxbin);
176   for(UShort_t det=1;det<=3;det++) {
177     //TObjArray* detInputArray = (TObjArray*)fInputArray->At(det);
178     TObjArray* detArray = (TObjArray*)fArray.At(det);
179     Int_t nRings = (det==1 ? 1 : 2);
180     for (UShort_t ir = 0; ir < nRings; ir++) {
181       Char_t ringChar = (ir == 0 ? 'I' : 'O');
182       //TObjArray* vtxInputArray = (TObjArray*)detInputArray->At(ir);
183       TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
184       TH2F* hMultTotal = (TH2F*)vtxArray->At(vtxbin);
185       
186       
187       TH2F* hMultInput = (TH2F*)fInputList->FindObject(Form("mult_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
188       
189       hMultTotal->Add(hMultInput);
190       
191       
192     }
193   }
194   
195   if(fMCevent && fPrimary)
196     ProcessPrimary();
197   
198   if(fStandalone) {
199     PostData(0, fOutputList); 
200   }
201   
202 }
203 //_____________________________________________________________________
204 void AliFMDAnalysisTaskDndeta::Terminate(Option_t */*option*/) {
205   
206   
207   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
208   
209   Int_t nVtxbins = pars->GetNvtxBins();
210   
211   for(UShort_t det=1;det<=3;det++) {
212     TObjArray* detArray = (TObjArray*)fArray.At(det);
213     Int_t nRings = (det==1 ? 1 : 2);
214     for (UShort_t ir = 0; ir < nRings; ir++) {
215       TObjArray* vtxArray = (TObjArray*)detArray->At(ir);
216       Char_t ringChar = (ir == 0 ? 'I' : 'O');
217       for(Int_t i =0; i<nVtxbins; i++) {
218         TH2F* hMultTotal = (TH2F*)vtxArray->At(i);
219         TH1D* hMultProj   = hMultTotal->ProjectionX(Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ringChar,i),1,hMultTotal->GetNbinsY());
220         fOutputList->Add(hMultProj);
221       }
222     }
223   }
224   
225 }
226 //_____________________________________________________________________
227 void AliFMDAnalysisTaskDndeta::ProcessPrimary() {
228   
229   fLastTrackByStrip.Reset(-1);
230   
231   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
232   
233   AliMCParticle* particle = 0;
234   AliStack* stack = fMCevent->Stack();
235   
236   TH1F* hPrimary = (TH1F*)fOutputList->FindObject("hMultvsEta");
237   AliHeader* header            = fMCevent->Header();
238   AliGenEventHeader* genHeader = header->GenEventHeader();
239   
240   TArrayF vertex;
241   genHeader->PrimaryVertex(vertex);
242   if(TMath::Abs(vertex.At(2)) > pars->GetVtxCutZ())
243     return;
244   Double_t delta           = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
245   Double_t vertexBinDouble = (vertex.At(2) + pars->GetVtxCutZ()) / delta;
246   Int_t    vertexBin       = (Int_t)vertexBinDouble;
247   
248   Bool_t firstTrack = kTRUE;
249   
250   // we loop over the primaries only unless we need the hits (diagnostics running slowly)
251   Int_t nTracks = stack->GetNprimary();
252   if(fRecordHits)
253     nTracks = stack->GetNtrack();
254   
255   for(Int_t i = 0 ;i<nTracks;i++) {
256     particle = fMCevent->GetTrack(i);
257     if(!particle)
258       continue;
259    
260     if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
261       hPrimary->Fill(particle->Eta());
262       
263
264       TH1F* hPrimVtxBin = (TH1F*)fOutputList->FindObject(Form("primmult_vtxbin%d",vertexBin));
265       hPrimVtxBin->Fill(particle->Eta());
266       if(firstTrack) {
267         fNMCevents.Fill(vertexBin);
268         firstTrack = kFALSE;
269       }
270     
271     }
272     if(fRecordHits) {
273       for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) {
274       
275       AliTrackReference* ref = particle->GetTrackReference(j);
276       UShort_t det,sec,strip;
277       Char_t   ring;
278       if(ref->DetectorId() != AliTrackReference::kFMD)
279         continue;
280       AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip);
281       Float_t thisStripTrack = fLastTrackByStrip(det,ring,sec,strip);
282       if(particle->Charge() != 0 && i != thisStripTrack ) {
283         //Double_t x,y,z;
284         
285         Float_t   eta   = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));//-1*TMath::Log(TMath::Tan(0.5*theta));
286         TH1F* hHits = (TH1F*)fOutputList->FindObject(Form("hHits_FMD%d%c_vtxbin%d",det,ring,vertexBin));
287         hHits->Fill(eta);
288         Float_t nstrips = (ring =='O' ? 256 : 512);
289         
290         fLastTrackByStrip(det,ring,sec,strip) = (Float_t)i;
291         
292         if(strip >0)
293           fLastTrackByStrip(det,ring,sec,strip-1) = (Float_t)i;
294         if(strip < (nstrips - 1))
295           fLastTrackByStrip(det,ring,sec,strip+1) = (Float_t)i;
296         
297         
298       }
299       }
300     }
301     
302   }
303   
304 }
305 //_____________________________________________________________________
306 //
307 //
308 // EOF