473b893b951bffe1f50d90d5a875105f9bc0b0be
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / 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 "AliGenPythiaEventHeader.h"
26 #include "AliHeader.h"
27 //#include "TDatabasePDG.h"
28 //#include "TParticlePDG.h"
29 #include "AliFMDStripIndex.h"
30 #include "AliESDInputHandler.h"
31
32 ClassImp(AliFMDAnalysisTaskDndeta)
33
34
35 AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta()
36 : fDebug(0),
37   fOutputList(0),
38   fInputList(0),
39   fVertexString(0x0),
40   fNevents(),
41   fNNSDevents(),
42   fNMCevents(),
43   fNMCNSDevents(),
44   fStandalone(kTRUE),
45   fLastTrackByStrip(0),
46   fVtxEff(1),
47   fVtxEffNSD(1)
48 {
49   // Default constructor
50   DefineInput (0, TList::Class());
51   DefineOutput(0, TList::Class());
52 }
53 //_____________________________________________________________________
54 AliFMDAnalysisTaskDndeta::AliFMDAnalysisTaskDndeta(const char* name, Bool_t SE):
55     AliAnalysisTask(name, "Density"),
56     fDebug(0),
57     fOutputList(0),
58     fInputList(0),
59     fVertexString(0x0),
60     fNevents(),
61     fNNSDevents(),
62     fNMCevents(),
63     fNMCNSDevents(),
64     fStandalone(kTRUE),
65     fLastTrackByStrip(0),
66     fVtxEff(1),
67     fVtxEffNSD(1)
68 {
69   fStandalone = SE;
70   if(fStandalone) {
71     DefineInput (0, TList::Class());
72     DefineInput(1, TObjString::Class());
73     DefineOutput(0, TList::Class());
74     
75   }
76 }
77 //_____________________________________________________________________
78 void AliFMDAnalysisTaskDndeta::CreateOutputObjects()
79 {
80   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
81   
82   if(!fOutputList)
83     fOutputList = new TList();
84   fOutputList->SetName("BackgroundCorrected");
85   
86   
87   TH2F* hMult = 0;
88   TH2F* hMultNSD = 0;
89   TH1F* hHits = 0;
90   TH2F* hMultTrVtx = 0;
91   TH1F* hPrimVertexBin = 0;
92   TH1F* hPrimVertexBinNSD = 0;
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   TH1F* hPrimaryNSD = new TH1F("hMultvsEtaNSD","hMultvsEtaNSD",
102                             hBgTmp->GetNbinsX(),
103                             hBgTmp->GetXaxis()->GetXmin(),
104                             hBgTmp->GetXaxis()->GetXmax());
105   hPrimaryNSD->Sumw2();
106   fOutputList->Add(hPrimaryNSD);
107   Int_t nVtxbins = pars->GetNvtxBins();
108   TH2F* hBg = 0;
109   for(Int_t i = 0; i< nVtxbins; i++) {
110      
111     for(Int_t det =1; det<=3;det++)
112       {
113         Int_t nRings = (det==1 ? 1 : 2);
114         for(Int_t ring = 0;ring<nRings;ring++)
115           {
116             Char_t ringChar = (ring == 0 ? 'I' : 'O');
117             Int_t  nSec     = (ring == 0 ? 20 : 40);
118             
119             
120             
121             hBg = pars->GetBackgroundCorrection(det, ringChar, i);
122             hMult  = new TH2F(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i),
123                               hBg->GetNbinsX(),
124                               hBg->GetXaxis()->GetXmin(),
125                               hBg->GetXaxis()->GetXmax(),
126                               nSec, 0, 2*TMath::Pi());
127             hMultTrVtx  = new TH2F(Form("dNdetaTrVtx_FMD%d%c_vtxbin%d",det,ringChar,i),Form("dNdetaTrVtx_FMD%d%c_vtxbin%d",det,ringChar,i),
128                                    hBg->GetNbinsX(),
129                                    hBg->GetXaxis()->GetXmin(),
130                                    hBg->GetXaxis()->GetXmax(),
131                                    nSec, 0, 2*TMath::Pi());
132             hMultNSD  = new TH2F(Form("dNdetaNSD_FMD%d%c_vtxbin%d",det,ringChar,i),Form("dNdetaNSD_FMD%d%c_vtxbin%d",det,ringChar,i),
133                                    hBg->GetNbinsX(),
134                                    hBg->GetXaxis()->GetXmin(),
135                                    hBg->GetXaxis()->GetXmax(),
136                                    nSec, 0, 2*TMath::Pi());
137             hHits  = new TH1F(Form("hMCHits_FMD%d%c_vtxbin%d",det,ringChar,i),Form("hMCHits_FMD%d%c_vtxbin%d",det,ringChar,i),
138                               hBg->GetNbinsX(),
139                               hBg->GetXaxis()->GetXmin(),
140                               hBg->GetXaxis()->GetXmax());
141             hHits->Sumw2();
142             fOutputList->Add(hHits);
143             
144             hMult->Sumw2();
145             fOutputList->Add(hMult);
146
147             hMultTrVtx->Sumw2();
148             fOutputList->Add(hMultTrVtx);
149             
150             hMultNSD->Sumw2();
151             fOutputList->Add(hMultNSD);
152             
153           }
154       } 
155   }
156   
157   for(Int_t i = 0; i< nVtxbins; i++) {
158    
159     hPrimVertexBin = new TH1F(Form("primmult_vtxbin%d",i),
160                               Form("primmult_vtxbin%d",i),
161                               hBgTmp->GetNbinsX(),
162                               hBgTmp->GetXaxis()->GetXmin(),
163                               hBgTmp->GetXaxis()->GetXmax());
164     hPrimVertexBin->Sumw2();
165     fOutputList->Add(hPrimVertexBin);
166     
167     hPrimVertexBinNSD = new TH1F(Form("primmult_NSD_vtxbin%d",i),
168                                  Form("primmult_NSD_vtxbin%d",i),
169                                  hBgTmp->GetNbinsX(),
170                                  hBgTmp->GetXaxis()->GetXmin(),
171                                  hBgTmp->GetXaxis()->GetXmax());
172     hPrimVertexBinNSD->Sumw2();
173     fOutputList->Add(hPrimVertexBinNSD);
174     
175
176     //SPD part
177     TH2F* hSPDMult  = new TH2F(Form("dNdeta_SPD_vtxbin%d",i),Form("dNdeta_SPD_vtxbin%d",i),
178                          hBgTmp->GetNbinsX(),
179                          hBgTmp->GetXaxis()->GetXmin(),
180                          hBgTmp->GetXaxis()->GetXmax(),
181                          20, 0, 2*TMath::Pi());
182     hSPDMult->Sumw2();
183     fOutputList->Add(hSPDMult);
184     TH2F* hSPDMultTrVtx  = new TH2F(Form("dNdetaTrVtx_SPD_vtxbin%d",i),Form("dNdetaTrVtx_SPD_vtxbin%d",i),
185                          hBgTmp->GetNbinsX(),
186                          hBgTmp->GetXaxis()->GetXmin(),
187                          hBgTmp->GetXaxis()->GetXmax(),
188                          20, 0, 2*TMath::Pi());
189     hSPDMultTrVtx->Sumw2();
190     fOutputList->Add(hSPDMultTrVtx);
191
192     TH2F* hSPDMultNSD  = new TH2F(Form("dNdetaNSD_SPD_vtxbin%d",i),Form("dNdetaNSD_SPD_vtxbin%d",i),
193                          hBgTmp->GetNbinsX(),
194                          hBgTmp->GetXaxis()->GetXmin(),
195                          hBgTmp->GetXaxis()->GetXmax(),
196                          20, 0, 2*TMath::Pi());
197     hSPDMultNSD->Sumw2();
198     fOutputList->Add(hSPDMultNSD);
199
200   }
201   
202   //AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
203   
204   //  TH2F* dNdetadphiHistogramTotal = new TH2F("dNdetadphiHistogramTotal","dNdetadphiHistogram;#eta;#Phi",pars->GetNetaBins(),-6,6,20,0,2*TMath::Pi());
205   //dNdetadphiHistogramTotal->SetErrorOption("g");
206   //fOutputList->Add(dNdetadphiHistogramTotal);
207   
208   
209   
210   fNevents.SetBins(nVtxbins,0,nVtxbins);
211   fNevents.SetName("nEvents");
212   fNNSDevents.SetBins(nVtxbins,0,nVtxbins);
213   fNNSDevents.SetName("nNSDEvents");
214   
215   fNMCevents.SetBins(nVtxbins,0,nVtxbins);
216   fNMCevents.SetName("nMCEvents");
217   fNMCNSDevents.SetBins(nVtxbins,0,nVtxbins);
218   fNMCNSDevents.SetName("nMCNSDEvents");
219   fOutputList->Add(&fNevents);
220   fOutputList->Add(&fNNSDevents);
221   fOutputList->Add(&fNMCevents);
222   fOutputList->Add(&fNMCNSDevents);
223 }
224 //_____________________________________________________________________
225 void AliFMDAnalysisTaskDndeta::ConnectInputData(Option_t */*option*/)
226 {
227   if(fStandalone) {
228     fInputList   = (TList*)GetInputData(0);
229   }
230 }
231 //_____________________________________________________________________
232 void AliFMDAnalysisTaskDndeta::Exec(Option_t */*option*/)
233 {
234   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
235   fVertexString = (TObjString*)fInputList->At(0);
236   
237   Int_t vtxbin   = fVertexString->GetString().Atoi();
238   
239   fNevents.Fill(vtxbin);
240   
241   //AliESDInputHandler* eventHandler = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
242   //AliESDEvent* esd = eventHandler->GetEvent();
243   Bool_t nsd = pars->IsEventTriggered(AliFMDAnaParameters::kNSD);
244   if(nsd) fNNSDevents.Fill(vtxbin);
245   
246   for(UShort_t det=1;det<=3;det++) {
247     Int_t nRings = (det==1 ? 1 : 2);
248     for (UShort_t ir = 0; ir < nRings; ir++) {
249       Char_t ringChar = (ir == 0 ? 'I' : 'O');
250             
251       TH2F* hMultTotal      = (TH2F*)fOutputList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
252       TH2F* hMultTotalTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
253       TH2F* hMultTotalNSD   = (TH2F*)fOutputList->FindObject(Form("dNdetaNSD_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
254       
255       
256       TH2F* hMultInput      = (TH2F*)fInputList->FindObject(Form("mult_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
257       TH2F* hMultInputTrVtx = (TH2F*)fInputList->FindObject(Form("multTrVtx_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
258       TH2F* hMultInputNSD   = (TH2F*)fInputList->FindObject(Form("multNSD_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
259       
260       hMultTotal->Add(hMultInput);
261       hMultTotalTrVtx->Add(hMultInputTrVtx);
262       if(nsd)
263         hMultTotalNSD->Add(hMultInputNSD);
264       
265     }
266   }
267   
268   //SPD code
269   TH2F* hMultSPDTotal      = (TH2F*)fOutputList->FindObject(Form("dNdeta_SPD_vtxbin%d",vtxbin));
270   TH2F* hMultSPDTotalTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_SPD_vtxbin%d",vtxbin));
271   TH2F* hMultSPDTotalNSD   = (TH2F*)fOutputList->FindObject(Form("dNdetaNSD_SPD_vtxbin%d",vtxbin));
272   TH2F* hMultSPDInput      = (TH2F*)fInputList->FindObject(Form("mult_SPD_vtxbin%d",vtxbin));
273   TH2F* hMultSPDInputTrVtx = (TH2F*)fInputList->FindObject(Form("multTrVtx_SPD_vtxbin%d",vtxbin));
274   TH2F* hMultSPDInputNSD   = (TH2F*)fInputList->FindObject(Form("multNSD_SPD_vtxbin%d",vtxbin));
275   
276   hMultSPDTotal->Add(hMultSPDInput);
277   hMultSPDTotalTrVtx->Add(hMultSPDInputTrVtx);
278   if(nsd)
279     hMultSPDTotalNSD->Add(hMultSPDInputNSD);
280   
281   if(pars->GetProcessPrimary())
282     ProcessPrimary();
283   
284   //TH2F* dNdetadphiHistogram = (TH2F*)fOutputList->FindObject("dNdetadphiHistogramSPDTrVtx");
285   
286   //  TH2F* dNdetadphiHistogramTotal = (TH2F*)fOutputList->FindObject("dNdetadphiHistogramTotal");
287   
288   //  if(vtxbin == 4)
289   //  dNdetadphiHistogramTotal->Add(dNdetadphiHistogram);
290   
291   
292   
293   
294   if(fStandalone) {
295     PostData(0, fOutputList); 
296   }
297   
298 }
299 //_____________________________________________________________________
300 void AliFMDAnalysisTaskDndeta::Terminate(Option_t */*option*/) {
301   
302   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
303   
304   Int_t nVtxbins = pars->GetNvtxBins();
305   
306   for(UShort_t det=1;det<=3;det++) {
307     Int_t nRings = (det==1 ? 1 : 2);
308     for (UShort_t ir = 0; ir < nRings; ir++) {
309       Char_t ringChar = (ir == 0 ? 'I' : 'O');
310       for(Int_t i =0; i<nVtxbins; i++) {
311         
312         TH2F* hMultTotal = (TH2F*)fOutputList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,i));
313         if(fVtxEff)
314           hMultTotal->Scale(fVtxEff);
315         
316         TH2F* hMultTotalNSD = (TH2F*)fOutputList->FindObject(Form("dNdetaNSD_FMD%d%c_vtxbin%d",det,ringChar,i));
317         if(fVtxEffNSD)
318           hMultTotalNSD->Scale(fVtxEffNSD);
319         
320         //TH2F* hMultTrVtx = (TH2F*)hMultTotal->Clone(Form("dNdeta_FMD%d%c_TrVtx_vtxbin%d",det,ringChar,i));
321         TH2F* hMultTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_FMD%d%c_vtxbin%d",det,ringChar,i));
322         
323         TH1D* hMultProj   = hMultTotal->ProjectionX(Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ringChar,i),1,hMultTotal->GetNbinsY());
324         TH1D* hMultProjTrVtx   = hMultTrVtx->ProjectionX(Form("dNdeta_FMD%d%c_TrVtx_vtxbin%d_proj",det,ringChar,i),1,hMultTrVtx->GetNbinsY());
325         TH1D* hMultProjNSD   = hMultTotalNSD->ProjectionX(Form("dNdetaNSD_FMD%d%c_vtxbin%d_proj",det,ringChar,i),1,hMultTotalNSD->GetNbinsY());
326         
327         //fOutputList->Add(hMultTrVtx);
328         fOutputList->Add(hMultProj);
329         fOutputList->Add(hMultProjTrVtx);
330         fOutputList->Add(hMultProjNSD);
331       }
332     }
333   }
334   
335   for(Int_t i =0; i<nVtxbins; i++) {
336     
337     TH2F* hSPDMult = (TH2F*)fOutputList->FindObject(Form("dNdeta_SPD_vtxbin%d",i));
338     TH2F* hSPDMultTrVtx = (TH2F*)fOutputList->FindObject(Form("dNdetaTrVtx_SPD_vtxbin%d",i));
339     TH2F* hSPDMultNSD = (TH2F*)fOutputList->FindObject(Form("dNdetaNSD_SPD_vtxbin%d",i));
340     
341     if(fVtxEff)
342       hSPDMult->Scale(fVtxEff);
343         
344     if(fVtxEffNSD)
345       hSPDMult->Scale(fVtxEffNSD);
346     
347     TH1D* hMultProj   = hSPDMult->ProjectionX(Form("dNdeta_SPD_vtxbin%d_proj",i),1,hSPDMult->GetNbinsY());
348     TH1D* hMultProjTrVtx   = hSPDMultTrVtx->ProjectionX(Form("dNdetaTrVtx_SPD_vtxbin%d_proj",i),1,hSPDMultTrVtx->GetNbinsY());
349     TH1D* hMultProjNSD   = hSPDMultNSD->ProjectionX(Form("dNdetaNSD_SPD_vtxbin%d_proj",i),1,hSPDMultNSD->GetNbinsY());
350     fOutputList->Add(hMultProj);
351     fOutputList->Add(hMultProjTrVtx);
352     fOutputList->Add(hMultProjNSD);
353   
354   }
355   
356   std::cout<<"FMD analysis accepted "<<fNevents.GetEntries()<<" INEL events and "<<fNNSDevents.GetEntries()<<" NSD events"<<std::endl;
357 }
358 //_____________________________________________________________________
359 void AliFMDAnalysisTaskDndeta::ProcessPrimary() {
360   
361   AliMCEventHandler* eventHandler = dynamic_cast<AliMCEventHandler*> (AliAnalysisManager::GetAnalysisManager()->GetMCtruthEventHandler());
362   AliMCEvent* mcEvent = eventHandler->MCEvent();
363   if(!mcEvent)
364     return;
365   
366   fLastTrackByStrip.Reset(-1);
367   
368   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
369   
370   AliMCParticle* particle = 0;
371   AliStack* stack = mcEvent->Stack();
372   
373   TH1F* hPrimary = (TH1F*)fOutputList->FindObject("hMultvsEta");
374   TH1F* hPrimaryNSD = (TH1F*)fOutputList->FindObject("hMultvsEtaNSD");
375   AliHeader* header            = mcEvent->Header();
376   AliGenEventHeader* genHeader = header->GenEventHeader();
377   
378   AliGenPythiaEventHeader* pythiaGenHeader = dynamic_cast<AliGenPythiaEventHeader*>(genHeader);
379   Bool_t nsd = kTRUE;
380   if (!pythiaGenHeader) {
381     std::cout<<" no pythia header!"<<std::endl;
382     //  return; 
383   }
384   else {
385     Int_t pythiaType = pythiaGenHeader->ProcessType();
386     
387     if(pythiaType==92||pythiaType==93)
388       nsd = kFALSE;
389   }
390   TArrayF vertex;
391   genHeader->PrimaryVertex(vertex);
392   if(TMath::Abs(vertex.At(2)) > pars->GetVtxCutZ())
393     return;
394   Double_t delta           = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
395   Double_t vertexBinDouble = (vertex.At(2) + pars->GetVtxCutZ()) / delta;
396   Int_t    vertexBin       = (Int_t)vertexBinDouble;
397     
398   Bool_t firstTrack = kTRUE;
399   Bool_t firstTrackNSD = kTRUE;
400   
401   TH1F* hPrimVtxBinNSD = (TH1F*)fOutputList->FindObject(Form("primmult_NSD_vtxbin%d",vertexBin));
402   TH1F* hPrimVtxBin = (TH1F*)fOutputList->FindObject(Form("primmult_vtxbin%d",vertexBin));
403   
404   // we loop over the primaries only unless we need the hits (diagnostics running slowly)
405   Int_t nTracks = stack->GetNprimary();
406   if(pars->GetProcessHits())
407     nTracks = stack->GetNtrack();
408   
409   for(Int_t i = 0 ;i<nTracks;i++) {
410     particle = (AliMCParticle*) mcEvent->GetTrack(i);
411     if(!particle)
412       continue;
413    
414     if(stack->IsPhysicalPrimary(i) && particle->Charge() != 0) {
415       hPrimary->Fill(particle->Eta());
416       hPrimVtxBin->Fill(particle->Eta());
417       if(firstTrack) {
418         fNMCevents.Fill(vertexBin);
419         firstTrack = kFALSE;
420       }
421       
422       if(nsd) {
423         hPrimaryNSD->Fill(particle->Eta());
424         hPrimVtxBinNSD->Fill(particle->Eta());
425         if(firstTrackNSD) {
426           fNMCNSDevents.Fill(vertexBin);
427           firstTrackNSD = kFALSE;
428         }
429       }
430       
431       
432     
433     }
434     if(pars->GetProcessHits()) {
435            
436       for(Int_t j=0; j<particle->GetNumberOfTrackReferences();j++) {
437         
438         AliTrackReference* ref = particle->GetTrackReference(j);
439         UShort_t det,sec,strip;
440         Char_t   ring;
441         if(ref->DetectorId() != AliTrackReference::kFMD)
442           continue;
443         AliFMDStripIndex::Unpack(ref->UserId(),det,ring,sec,strip);
444         Float_t thisStripTrack = fLastTrackByStrip(det,ring,sec,strip);
445         if(particle->Charge() != 0 && i != thisStripTrack ) {
446           //Double_t x,y,z;
447           
448           Float_t   eta   = pars->GetEtaFromStrip(det,ring,sec,strip,vertex.At(2));//-1*TMath::Log(TMath::Tan(0.5*theta));
449           TH1F* hHits = (TH1F*)fOutputList->FindObject(Form("hMCHits_FMD%d%c_vtxbin%d",det,ring,vertexBin));
450           
451         
452           hHits->Fill(eta);
453           
454           Float_t nstrips = (ring =='O' ? 256 : 512);
455           
456           fLastTrackByStrip(det,ring,sec,strip) = (Float_t)i;
457         
458           if(strip >0)
459             fLastTrackByStrip(det,ring,sec,strip-1) = (Float_t)i;
460           if(strip < (nstrips - 1))
461             fLastTrackByStrip(det,ring,sec,strip+1) = (Float_t)i;
462           
463         }
464       
465         
466       }
467       
468       
469     }
470
471       
472       
473   }
474 }
475 //_____________________________________________________________________
476 //
477 //
478 // EOF