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