All FMD corrections are now divided into the analysis + adding new corrections
[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   fVertexString(0x0),
37   fNevents(),
38   fNMCevents(),
39   fStandalone(kTRUE),
40   fLastTrackByStrip(0)
41 {
42   // Default constructor
43   DefineInput (0, TList::Class());
44   DefineOutput(0, TList::Class());
45 }
46 //_____________________________________________________________________
47 AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta(const char* name, Bool_t SE):
48     AliAnalysisTask(name, "Density"),
49     fDebug(0),
50     fOutputList(0),
51     fInputList(0),
52     fVertexString(0x0),
53     fNevents(),
54     fNMCevents(),
55     fStandalone(kTRUE),
56     fLastTrackByStrip(0)
57 {
58   fStandalone = SE;
59   if(fStandalone) {
60     DefineInput (0, TList::Class());
61     DefineInput(1, TObjString::Class());
62     DefineOutput(0, TList::Class());
63     
64   }
65 }
66 //_____________________________________________________________________
67 void AliFMDAnalysisTaskDndeta::CreateOutputObjects()
68 {
69   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
70   
71   if(!fOutputList)
72     fOutputList = new TList();
73   fOutputList->SetName("BackgroundCorrected");
74   
75   
76   TH2F* hMult = 0;
77   TH1F* hHits = 0;
78   TH1F* hPrimVertexBin = 0;
79   
80   
81   TH2F* hBgTmp   = pars->GetBackgroundCorrection(1, 'I', 0);
82   TH1F* hPrimary = new TH1F("hMultvsEta","hMultvsEta",
83                             hBgTmp->GetNbinsX(),
84                             hBgTmp->GetXaxis()->GetXmin(),
85                             hBgTmp->GetXaxis()->GetXmax());
86   hPrimary->Sumw2();
87   fOutputList->Add(hPrimary);
88   Int_t nVtxbins = pars->GetNvtxBins();
89   TH2F* hBg = 0;
90   for(Int_t i = 0; i< nVtxbins; i++) {
91      
92     for(Int_t det =1; det<=3;det++)
93       {
94         Int_t nRings = (det==1 ? 1 : 2);
95         for(Int_t ring = 0;ring<nRings;ring++)
96           {
97             Char_t ringChar = (ring == 0 ? 'I' : 'O');
98             Int_t  nSec     = (ring == 0 ? 20 : 40);
99             
100             
101             
102             hBg = pars->GetBackgroundCorrection(det, ringChar, i);
103             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->GetNbinsX(),
105                               hBg->GetXaxis()->GetXmin(),
106                               hBg->GetXaxis()->GetXmax(),
107                               nSec, 0, 2*TMath::Pi());
108             
109             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->GetNbinsX(),
111                               hBg->GetXaxis()->GetXmin(),
112                               hBg->GetXaxis()->GetXmax());
113             hHits->Sumw2();
114             fOutputList->Add(hHits);
115             
116             hMult->Sumw2();
117             fOutputList->Add(hMult);
118             
119           }
120       } 
121   }
122   
123   for(Int_t i = 0; i< nVtxbins; i++) {
124    
125     hPrimVertexBin = new TH1F(Form("primmult_vtxbin%d",i),
126                               Form("primmult_vtxbin%d",i),
127                               hBgTmp->GetNbinsX(),
128                               hBgTmp->GetXaxis()->GetXmin(),
129                               hBgTmp->GetXaxis()->GetXmax());
130     hPrimVertexBin->Sumw2();
131     fOutputList->Add(hPrimVertexBin);
132     
133   }
134   
135   fNevents.SetBins(nVtxbins,0,nVtxbins);
136   fNevents.SetName("nEvents");
137   fNMCevents.SetBins(nVtxbins,0,nVtxbins);
138   fNMCevents.SetName("nMCEvents");
139   fOutputList->Add(&fNevents);
140   fOutputList->Add(&fNMCevents);
141   
142 }
143 //_____________________________________________________________________
144 void AliFMDAnalysisTaskDndeta::ConnectInputData(Option_t */*option*/)
145 {
146   if(fStandalone) {
147     fInputList   = (TList*)GetInputData(0);
148   }
149 }
150 //_____________________________________________________________________
151 void AliFMDAnalysisTaskDndeta::Exec(Option_t */*option*/)
152 {
153   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
154   fVertexString = (TObjString*)fInputList->At(0);
155   
156   Int_t vtxbin   = fVertexString->GetString().Atoi();
157   
158   fNevents.Fill(vtxbin);
159   
160   for(UShort_t det=1;det<=3;det++) {
161     Int_t nRings = (det==1 ? 1 : 2);
162     for (UShort_t ir = 0; ir < nRings; ir++) {
163       Char_t ringChar = (ir == 0 ? 'I' : 'O');
164             
165       TH2F* hMultTotal = (TH2F*)fOutputList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
166      
167       
168       TH2F* hMultInput = (TH2F*)fInputList->FindObject(Form("mult_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
169       
170       hMultTotal->Add(hMultInput);
171             
172     }
173   }
174   
175   if(pars->GetProcessPrimary())
176     ProcessPrimary();
177   
178   if(fStandalone) {
179     PostData(0, fOutputList); 
180   }
181   
182 }
183 //_____________________________________________________________________
184 void AliFMDAnalysisTaskDndeta::Terminate(Option_t */*option*/) {
185   
186   
187   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
188   
189   Int_t nVtxbins = pars->GetNvtxBins();
190   
191   for(UShort_t det=1;det<=3;det++) {
192     Int_t nRings = (det==1 ? 1 : 2);
193     for (UShort_t ir = 0; ir < nRings; ir++) {
194       Char_t ringChar = (ir == 0 ? 'I' : 'O');
195       for(Int_t i =0; i<nVtxbins; i++) {
196         
197         TH2F* hMultTotal = (TH2F*)fOutputList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i));
198         TH2F* hMultTrVtx = (TH2F*)hMultTotal->Clone(Form("dNdeta_FMD%d%c_TrVtx_vtxbin%d",det,ringChar,i));
199         
200         hMultTotal->Scale(1/pars->GetEventSelectionEfficiency(i));
201                 
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);
207       }
208     }
209   }
210   
211 }
212 //_____________________________________________________________________
213 void AliFMDAnalysisTaskDndeta::ProcessPrimary() {
214   
215   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
216   AliMCEvent* mcEvent = eventHandler->MCEvent();
217   if(!mcEvent)
218     return;
219   
220   fLastTrackByStrip.Reset(-1);
221   
222   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
223   
224   AliMCParticle* particle = 0;
225   AliStack* stack = mcEvent->Stack();
226   
227   TH1F* hPrimary = (TH1F*)fOutputList->FindObject("hMultvsEta");
228   AliHeader* header            = mcEvent->Header();
229   AliGenEventHeader* genHeader = header->GenEventHeader();
230   
231   TArrayF vertex;
232   genHeader->PrimaryVertex(vertex);
233   if(TMath::Abs(vertex.At(2)) > pars->GetVtxCutZ())
234     return;
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;
238   
239   Bool_t firstTrack = kTRUE;
240   
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();
245   
246   for(Int_t i = 0 ;i<nTracks;i++) {
247     particle = (AliMCParticle*) mcEvent->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     if(pars->GetProcessHits()) {
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(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("hMCHits_FMD%d%c_vtxbin%d",det,ring,vertexBin));
279           
280         
281           hHits->Fill(eta);
282           
283           Float_t nstrips = (ring =='O' ? 256 : 512);
284           
285           fLastTrackByStrip(det,ring,sec,strip) = (Float_t)i;
286         
287           if(strip >0)
288             fLastTrackByStrip(det,ring,sec,strip-1) = (Float_t)i;
289           if(strip < (nstrips - 1))
290             fLastTrackByStrip(det,ring,sec,strip+1) = (Float_t)i;
291           
292         }
293       
294         
295       }
296       
297       
298     }
299
300       
301       
302   }
303 }
304 //_____________________________________________________________________
305 //
306 //
307 // EOF