Various upgrades and bug fixes etc
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis / AliFMDAnalysisTaskBackgroundCorrection.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 "TH2F.h"
10 #include "AliFMDAnalysisTaskBackgroundCorrection.h"
11 #include "AliAnalysisManager.h"
12 #include "AliESDFMD.h"
13 #include "AliESDEvent.h"
14 #include "AliAODEvent.h"
15 #include "AliAODHandler.h"
16 #include "AliMCEventHandler.h"
17 #include "AliStack.h"
18 #include "AliLog.h"
19 #include "AliESDVertex.h"
20 #include "TMath.h"
21 #include "AliFMDAnaParameters.h"
22 #include "AliESDInputHandler.h"
23 #include "AliMultiplicity.h"
24 #include "TProfile2D.h"
25 //#include "AliFMDGeometry.h"
26
27 ClassImp(AliFMDAnalysisTaskBackgroundCorrection)
28
29
30 AliFMDAnalysisTaskBackgroundCorrection::AliFMDAnalysisTaskBackgroundCorrection()
31 : fDebug(0),
32   fOutputList(0),
33   fInputList(0),
34   fHitList(0),
35   fVertexString(0x0),
36   fNevents(),
37   fStandalone(kTRUE),
38   fOutputVertexString(0)
39 {
40   // Default constructor
41   DefineInput (0, TList::Class());
42   DefineOutput(0, TList::Class());
43
44 }
45 //_____________________________________________________________________
46 AliFMDAnalysisTaskBackgroundCorrection::AliFMDAnalysisTaskBackgroundCorrection(const char* name, Bool_t SE):
47     AliAnalysisTask(name, "Density"),
48     fDebug(0),
49     fOutputList(0),
50     fInputList(0),
51     fHitList(0),
52     fVertexString(0x0),
53     fNevents(),
54     fStandalone(kTRUE),
55     fOutputVertexString(0)
56 {
57   fStandalone = SE;
58   if(fStandalone) {
59     DefineInput (0, TList::Class());
60     DefineOutput(0, TList::Class());
61   }
62 }
63 //_____________________________________________________________________
64 void AliFMDAnalysisTaskBackgroundCorrection::CreateOutputObjects()
65 {
66   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
67   
68   if(!fOutputList)
69     fOutputList = new TList();
70   fOutputList->SetName("BackgroundCorrectedPerEvent");
71   if(!fHitList)
72     fHitList = new TList();
73   fHitList->SetName("HitsList");
74   
75   //if(fStandalone) {
76   fOutputVertexString = new TObjString();
77   // }
78   fOutputList->Add(fOutputVertexString);
79   
80   
81   
82   TH2F* hMult         = 0;
83   TH2F* hMultTrVtx    = 0;
84   TH2F* hMultNSD      = 0;
85   TH2F* hHits         = 0;
86   TH2F* hSPDMult      = 0;
87   TH2F* hSPDMultTrVtx = 0;
88   TH2F* hSPDMultNSD   = 0;
89   // TH2F* hHitsNoCuts = 0;
90   Int_t nVtxbins = pars->GetNvtxBins();
91   for(Int_t i = 0; i< nVtxbins; i++) {
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             TH2F* hBg = pars->GetBackgroundCorrection(det, ringChar, i);
102             hMult  = new TH2F(Form("mult_FMD%d%c_vtxbin%d",det,ringChar,i),Form("mult_FMD%d%c_vtxbin%d",det,ringChar,i),
103                               hBg->GetNbinsX(),
104                               hBg->GetXaxis()->GetXmin(),
105                               hBg->GetXaxis()->GetXmax(),
106                               nSec, 0, 2*TMath::Pi());
107             hMult->Sumw2();
108             fOutputList->Add(hMult);
109             hMultTrVtx  = new TH2F(Form("multTrVtx_FMD%d%c_vtxbin%d",det,ringChar,i),Form("multTrVtx_FMD%d%c_vtxbin%d",det,ringChar,i),
110                                    hBg->GetNbinsX(),
111                                    hBg->GetXaxis()->GetXmin(),
112                                    hBg->GetXaxis()->GetXmax(),
113                                    nSec, 0, 2*TMath::Pi());
114             hMultTrVtx->Sumw2();
115
116             fOutputList->Add(hMultTrVtx);
117             
118             hMultNSD  = new TH2F(Form("multNSD_FMD%d%c_vtxbin%d",det,ringChar,i),Form("multNSD_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             hMultNSD->Sumw2();
124
125             fOutputList->Add(hMultNSD);
126             
127             hHits  = new TH2F(Form("hits_FMD%d%c_vtxbin%d",det,ringChar,i),Form("hits_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             
133             /*  hHitsNoCuts  = new TH2F(Form("hits_NoCuts_FMD%d%c_vtxbin%d",det,ringChar,i),Form("hits_NoCuts_FMD%d%c_vtxbin%d",det,ringChar,i),
134                                     hBg->GetNbinsX(),
135                                     hBg->GetXaxis()->GetXmin(),
136                                     hBg->GetXaxis()->GetXmax(),
137                                     nSec, 0, 2*TMath::Pi());
138             
139             */
140             hHits->Sumw2();
141             fHitList->Add(hHits);
142             
143           }
144       }
145     //HHD SPD hists
146     TH2F* hBg = pars->GetBackgroundCorrection(1, 'I', i);
147     hSPDMult  = new TH2F(Form("mult_SPD_vtxbin%d",i),Form("mult_SPD_vtxbin%d",i),
148                          hBg->GetNbinsX(),
149                          hBg->GetXaxis()->GetXmin(),
150                          hBg->GetXaxis()->GetXmax(),
151                          20, 0, 2*TMath::Pi());
152     hSPDMult->Sumw2();
153     fOutputList->Add(hSPDMult);
154     hSPDMultTrVtx  = new TH2F(Form("multTrVtx_SPD_vtxbin%d",i),Form("multTrVtx_SPD_vtxbin%d",i),
155                          hBg->GetNbinsX(),
156                          hBg->GetXaxis()->GetXmin(),
157                          hBg->GetXaxis()->GetXmax(),
158                          20, 0, 2*TMath::Pi());
159     hSPDMultTrVtx->Sumw2();
160     fOutputList->Add(hSPDMultTrVtx);
161     hSPDMultNSD  = new TH2F(Form("multNSD_SPD_vtxbin%d",i),Form("multNSD_SPD_vtxbin%d",i),
162                          hBg->GetNbinsX(),
163                          hBg->GetXaxis()->GetXmin(),
164                          hBg->GetXaxis()->GetXmax(),
165                          20, 0, 2*TMath::Pi());
166     hSPDMultNSD->Sumw2();
167     fOutputList->Add(hSPDMultNSD);
168     
169     
170   }
171   TH2F* hBg = pars->GetBackgroundCorrection(1, 'I', 5);
172   TH2F* dNdetadphiHistogram = new TH2F("dNdetadphiHistogramTrVtx","dNdetadphiHistogramTrVtx;#eta;#Phi",pars->GetNetaBins(),hBg->GetXaxis()->GetXmin(),hBg->GetXaxis()->GetXmax(),20,0,2*TMath::Pi());
173   TH2F* dNdetadphiHistogramSPD = new TH2F("dNdetadphiHistogramSPDTrVtx","dNdetadphiHistogramSPDTrVtx;#eta;#Phi",pars->GetNetaBins(),hBg->GetXaxis()->GetXmin(),hBg->GetXaxis()->GetXmax(),20,0,2*TMath::Pi());
174   //dNdetadphiHistogram->SetErrorOption("g");
175   
176   fHitList->Add(dNdetadphiHistogram);
177   fOutputList->Add(dNdetadphiHistogram);
178   fHitList->Add(dNdetadphiHistogramSPD);
179   fOutputList->Add(dNdetadphiHistogramSPD);
180   
181   
182 }
183 //_____________________________________________________________________
184 void AliFMDAnalysisTaskBackgroundCorrection::ConnectInputData(Option_t */*option*/)
185 {
186   if(fStandalone) {
187     fInputList   = (TList*)GetInputData(0);
188     
189   }
190 }
191 //_____________________________________________________________________
192 void AliFMDAnalysisTaskBackgroundCorrection::Exec(Option_t */*option*/)
193 {
194   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
195   
196   fVertexString = (TObjString*)fInputList->At(0);
197    
198   Int_t vtxbin   = fVertexString->GetString().Atoi();
199   fOutputVertexString->SetString(Form("%d",vtxbin));
200   
201   fNevents.Fill(vtxbin);
202   //Reset everything
203   for(UShort_t det=1;det<=3;det++) {
204     Int_t nRings = (det==1 ? 1 : 2);
205     for (UShort_t ir = 0; ir < nRings; ir++) {
206       Char_t ringChar = (ir == 0 ? 'I' : 'O');
207       TH2F* hMult = (TH2F*)fOutputList->FindObject(Form("mult_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
208       hMult->Reset();
209       TH2F* hMultTrVtx = (TH2F*)fOutputList->FindObject(Form("multTrVtx_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
210       hMultTrVtx->Reset();
211       TH2F* hMultNSD = (TH2F*)fOutputList->FindObject(Form("multNSD_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
212       hMultNSD->Reset();
213     
214       TH2F* hSPDMult      = (TH2F*)fOutputList->FindObject(Form("mult_SPD_vtxbin%d",vtxbin));
215       hSPDMult->Reset();
216       TH2F* hSPDMultTrVtx = (TH2F*)fOutputList->FindObject(Form("multTrVtx_SPD_vtxbin%d",vtxbin));
217       hSPDMultTrVtx->Reset();
218       TH2F* hSPDMultNSD = (TH2F*)fOutputList->FindObject(Form("multNSD_SPD_vtxbin%d",vtxbin));
219       hSPDMultNSD->Reset();
220     }
221     
222   }
223   
224  
225   Bool_t nsd = pars->IsEventTriggered(AliFMDAnaParameters::kNSD);
226   for(UShort_t det=1;det<=3;det++) {
227     
228     Int_t nRings = (det==1 ? 1 : 2);
229     for (UShort_t ir = 0; ir < nRings; ir++) {
230       Char_t ringChar = (ir == 0 ? 'I' : 'O');
231       
232       TH2F* hMult      = (TH2F*)fOutputList->FindObject(Form("mult_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
233       TH2F* hMultTrVtx = (TH2F*)fOutputList->FindObject(Form("multTrVtx_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
234       TH2F* hMultNSD   = (TH2F*)fOutputList->FindObject(Form("multNSD_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
235       TH2F* hMultInput = (TH2F*)fInputList->FindObject(Form("FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
236       TH2F* hHits      = (TH2F*)fHitList->FindObject(Form("hits_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
237       
238       //if(pars->GetProcessHits())
239       hHits->Add(hMultInput);
240       
241       TH2F* hBg        = pars->GetBackgroundCorrection(det, ringChar, vtxbin);
242       TH2F* hBgNSD     = pars->GetBackgroundCorrectionNSD(det, ringChar, vtxbin);
243       hMult->Add(hMultInput);
244       hMultTrVtx->Add(hMultInput);
245       
246       hMult->Divide(hBg);//,"B");
247       hMultTrVtx->Divide(hBg);//,"B");
248       
249       if(nsd) {
250         hMultNSD->Add(hMultInput);
251         hMultNSD->Divide(hBg);
252       }
253       
254       //sharing efficiency correction ?
255       if(pars->SharingEffPresent()) {
256         TH1F* hSharingEff = pars->GetSharingEfficiencyTrVtx(det,ringChar,vtxbin); 
257         TH1F* hSharingEffTrVtx = pars->GetSharingEfficiencyTrVtx(det,ringChar,vtxbin);  
258       
259         for(Int_t nx=1; nx<hMult->GetNbinsX(); nx++) {
260           Float_t correction      = hSharingEff->GetBinContent(nx);
261           Float_t correctionTrVtx = hSharingEffTrVtx->GetBinContent(nx);
262           //FIXME : This should be for NSD events
263           Float_t correctionNSD   = hSharingEff->GetBinContent(nx);
264           
265           for(Int_t ny=1; ny<hMult->GetNbinsY(); ny++) {
266             
267             if(correction != 0){
268               hMult->SetBinContent(nx,ny,hMult->GetBinContent(nx,ny)/correction);
269               Float_t error = TMath::Sqrt(TMath::Power(hMult->GetBinError(nx,ny),2) + TMath::Power(hMult->GetBinContent(nx,ny)*hSharingEff->GetBinError(nx),2)) / correction;
270               hMult->SetBinError(nx,ny,error);
271             }
272             if(correctionTrVtx != 0){
273               hMultTrVtx->SetBinContent(nx,ny,hMultTrVtx->GetBinContent(nx,ny)/correctionTrVtx);
274               Float_t error = TMath::Sqrt(TMath::Power(hMultTrVtx->GetBinError(nx,ny),2) + TMath::Power(hMultTrVtx->GetBinContent(nx,ny)*hSharingEffTrVtx->GetBinError(nx),2)) / correctionTrVtx;
275               hMultTrVtx->SetBinError(nx,ny,error);
276             }
277             if(correctionNSD != 0 && nsd) {
278               hMultNSD->SetBinContent(nx,ny,hMultNSD->GetBinContent(nx,ny)/correctionNSD);
279               Float_t error = TMath::Sqrt(TMath::Power(hMultNSD->GetBinError(nx,ny),2) + TMath::Power(hMultNSD->GetBinContent(nx,ny)*hSharingEff->GetBinError(nx),2)) / correctionNSD;
280               hMultNSD->SetBinError(nx,ny,error);
281             }
282             
283             
284           }
285           
286         }
287       }
288       
289       //if(pars->GetEventSelectionEfficiency(vtxbin) > 0)
290       //        hMult->Scale(1/pars->GetEventSelectionEfficiency(vtxbin));
291         //else
292         //hMult->Scale(0);
293       hMult->Divide(pars->GetEventSelectionEfficiency("INEL",vtxbin,ringChar));
294       hMultNSD->Divide(pars->GetEventSelectionEfficiency("NSD",vtxbin,ringChar));
295       
296       
297       }
298   }
299   
300   //HHD SPD code
301   
302   TH2F* hSPDMult      = (TH2F*)fOutputList->FindObject(Form("mult_SPD_vtxbin%d",vtxbin));
303   TH2F* hSPDMultTrVtx = (TH2F*)fOutputList->FindObject(Form("multTrVtx_SPD_vtxbin%d",vtxbin));
304   TH2F* hSPDMultNSD   = (TH2F*)fOutputList->FindObject(Form("multNSD_SPD_vtxbin%d",vtxbin));
305   
306   AliESDInputHandler* eventHandler = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
307   AliESDEvent* esd = eventHandler->GetEvent();
308   
309   const AliMultiplicity* spdmult = esd->GetMultiplicity();
310   for(Int_t j = 0; j< spdmult->GetNumberOfTracklets();j++) {
311     hSPDMult->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
312     hSPDMultTrVtx->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
313     hSPDMultNSD->Fill(spdmult->GetEta(j),spdmult->GetPhi(j));
314   }
315   for(Int_t j = 0; j< spdmult->GetNumberOfSingleClusters();j++) {
316     hSPDMult->Fill(-TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.)),spdmult->GetPhiSingle(j));
317     hSPDMultTrVtx->Fill(-TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.)),spdmult->GetPhiSingle(j));
318     hSPDMultNSD->Fill(-TMath::Log(TMath::Tan(spdmult->GetThetaSingle(j)/2.)),spdmult->GetPhiSingle(j));
319   }
320   
321   TH2F* hBgSPD        = pars->GetBackgroundCorrection(0, 'Q', vtxbin);
322   if(hBgSPD) { 
323    
324     TH1F* hDead      = pars->GetSPDDeadCorrection(vtxbin);
325     for(Int_t i=1; i<=hSPDMult->GetNbinsX(); i++) {
326       for(Int_t j=1; j<=hSPDMult->GetNbinsY(); j++) {
327         Float_t mult = hSPDMult->GetBinContent(i,j);
328         Float_t correction = hBgSPD->GetBinContent(i,j);
329         Float_t correctedMult = 0;
330         Float_t correctedError = 0;
331         if(correction > 0 && mult > 0) {
332           //correctedMult = mult;///correction;
333           correctedMult = mult/correction;
334           
335           if(hDead->GetBinContent(i) > 0)
336             correctedMult = correctedMult/hDead->GetBinContent(i);
337           correctedError = correctedMult*TMath::Sqrt( TMath::Power(hSPDMult->GetBinError(i,j)/hSPDMult->GetBinContent(i,j),2) + 
338                                                       TMath::Power(hBgSPD->GetBinError(i,j)/hBgSPD->GetBinContent(i,j),2)); 
339           
340             }
341         
342          if(correctedMult != 0) {
343            hSPDMult->SetBinContent(i,j,correctedMult);
344            hSPDMultTrVtx->SetBinContent(i,j,correctedMult);
345            hSPDMult->SetBinError(i,j,correctedError);
346            hSPDMultTrVtx->SetBinError(i,j,correctedError);
347            if(nsd) {
348              hSPDMultNSD->SetBinContent(i,j,correctedMult);
349              hSPDMultNSD->SetBinError(i,j,correctedError);
350         }
351            
352         }
353     }
354     }
355     hSPDMult->Divide(pars->GetEventSelectionEfficiency("INEL",vtxbin,'I'));
356     hSPDMultNSD->Divide(pars->GetEventSelectionEfficiency("NSD",vtxbin,'I'));
357   }
358   else
359     AliWarning("No SPD background map found");
360   
361   //std::cout<<spdmult->GetNumberOfTracklets()<<"  "<<spdmult->GetNumberOfITSClusters(0)<<"    "<< spdmult->GetNumberOfSingleClusters()<<std::endl;
362   
363   CreatePerEventHistogram(vtxbin);
364   
365   if(fStandalone) {
366     PostData(0, fOutputList); 
367   }
368   
369 }
370 //_____________________________________________________________________
371 void AliFMDAnalysisTaskBackgroundCorrection::Terminate(Option_t */*option*/) {
372   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
373   
374   Int_t nVtxbins = pars->GetNvtxBins();
375   
376   for(UShort_t det=1;det<=3;det++) {
377     Int_t nRings = (det==1 ? 1 : 2);
378     for (UShort_t ir = 0; ir < nRings; ir++) {
379       Char_t ringChar = (ir == 0 ? 'I' : 'O');
380       for(Int_t i =0; i<nVtxbins; i++) {
381         TH2F* hHits      = (TH2F*)fHitList->FindObject(Form("hits_FMD%d%c_vtxbin%d",det,ringChar,i));
382         TH1D* hHitsproj  = hHits->ProjectionX(Form("hits_FMD%d%c_vtxbin%d_proj",det,ringChar,i),1,hHits->GetNbinsY());
383         TH1D* hHitsNoCuts = (TH1D*)hHitsproj->Clone(Form("hits_NoCuts_FMD%d%c_vtxbin%d_proj",det,ringChar,i));
384         if(pars->GetEventSelectionEfficiency(i) > 0)
385           hHitsNoCuts->Scale(1/pars->GetEventSelectionEfficiency(i));
386         else
387           hHitsNoCuts->Scale(0);
388         fHitList->Add(hHitsproj);
389         fHitList->Add(hHitsNoCuts);
390         
391       }
392     }
393   }
394 }
395
396 //_____________________________________________________________________
397 void AliFMDAnalysisTaskBackgroundCorrection::CreatePerEventHistogram(Int_t vtxbin) {
398   
399   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
400   TH2F* dNdetadphiHistogram = (TH2F*)fHitList->FindObject("dNdetadphiHistogramTrVtx");
401   TH2F* dNdetadphiHistogramSPD = (TH2F*)fHitList->FindObject("dNdetadphiHistogramSPDTrVtx");
402   
403   dNdetadphiHistogram->Reset();
404   dNdetadphiHistogramSPD->Reset();
405   
406   for(Int_t det = 0; det<=3; det++) {
407     Int_t maxRing = (det<= 1 ? 0 : 1);
408     
409     
410     
411     for(Int_t ring = 0;ring<=maxRing;ring++) {
412       
413       Char_t ringChar = (ring == 0 ? 'I' : 'O');
414       
415       TH2F* multhistoriginal = 0;
416       
417       if(det == 0)
418         multhistoriginal = (TH2F*)fOutputList->FindObject(Form("multTrVtx_SPD_vtxbin%d",vtxbin));
419       else       
420         multhistoriginal = (TH2F*)fOutputList->FindObject(Form("multTrVtx_FMD%d%c_vtxbin%d",det,ringChar,vtxbin));
421       
422       TH2F* multhist = (TH2F*)multhistoriginal->Clone("tmp");
423       
424       
425       
426       if(ringChar == 'O' && det > 0)
427         multhist->RebinY(2);
428       
429       for(Int_t i=pars->GetFirstEtaBinToInclude(vtxbin,det,ringChar); i<=pars->GetLastEtaBinToInclude(vtxbin,det,ringChar); i++) {
430         for(Int_t j=1; j<=multhist->GetNbinsY(); j++) {
431           if(multhist->GetBinContent(i,j) < 0.0001) continue;
432           
433           Bool_t overlapFMD = kFALSE;
434           Bool_t overlapSPD = kFALSE;
435           
436           if(det == 1 && ringChar =='I')
437             if(i<=pars->GetLastEtaBinToInclude(vtxbin,2,'I'))
438               overlapFMD = kTRUE;
439                   
440           if(det == 2 && ringChar =='O') {
441             if(i>=pars->GetFirstEtaBinToInclude(vtxbin,2,'I'))
442               overlapFMD = kTRUE;
443             if(i<=pars->GetLastEtaBinToInclude(vtxbin,0,'I'))// && TMath::Abs(multhist->GetXaxis()->GetBinCenter(i)) < 2)
444               overlapSPD = kTRUE;
445           }
446           if(det == 2 && ringChar =='I')
447             if(i<=pars->GetLastEtaBinToInclude(vtxbin,2,'O') || i>=pars->GetFirstEtaBinToInclude(vtxbin,1,'I'))
448               overlapFMD = kTRUE;
449                   
450           if(det == 3 && ringChar =='I')
451             if(i>=pars->GetFirstEtaBinToInclude(vtxbin,3,'O'))
452               overlapFMD = kTRUE; 
453           
454           if(det == 3 && ringChar =='O') {
455             if(i<=pars->GetLastEtaBinToInclude(vtxbin,3,'I'))
456               overlapFMD = kTRUE;
457             if(i>=pars->GetFirstEtaBinToInclude(vtxbin,0,'I'))// && TMath::Abs(multhist->GetXaxis()->GetBinCenter(i)) < 2)
458               overlapSPD = kTRUE;
459           }
460           
461           if(det == 0) {
462             if(i<=pars->GetLastEtaBinToInclude(vtxbin,3,'O') || i>=pars->GetFirstEtaBinToInclude(vtxbin,2,'O'))
463               overlapSPD = kTRUE;
464           }
465           //std::cout<<overlapFMD<<"    "<<overlapSPD<<std::endl;
466           
467           
468           if(det > 0) {
469             if(overlapFMD) {
470               dNdetadphiHistogram->SetBinContent(i,j,dNdetadphiHistogram->GetBinContent(i,j)+0.5*multhist->GetBinContent(i,j));
471               dNdetadphiHistogram->SetBinError(i,j,TMath::Sqrt(TMath::Power(dNdetadphiHistogram->GetBinError(i,j),2)+TMath::Power(0.5*multhist->GetBinError(i,j),2)));
472             }
473             else {
474               dNdetadphiHistogram->SetBinContent(i,j,multhist->GetBinContent(i,j));
475               dNdetadphiHistogram->SetBinError(i,j,multhist->GetBinError(i,j));
476             }
477           }
478           
479           if( overlapFMD && overlapSPD) {
480             dNdetadphiHistogramSPD->SetBinContent(i,j,dNdetadphiHistogramSPD->GetBinContent(i,j)+0.33*multhist->GetBinContent(i,j));
481             dNdetadphiHistogramSPD->SetBinError(i,j,TMath::Sqrt(TMath::Power(dNdetadphiHistogramSPD->GetBinError(i,j),2)+TMath::Power(0.33*multhist->GetBinError(i,j),2)));
482           }
483           else if( overlapFMD || overlapSPD) {
484             dNdetadphiHistogramSPD->SetBinContent(i,j,dNdetadphiHistogramSPD->GetBinContent(i,j)+0.5*multhist->GetBinContent(i,j));
485             dNdetadphiHistogramSPD->SetBinError(i,j,TMath::Sqrt(TMath::Power(dNdetadphiHistogramSPD->GetBinError(i,j),2)+TMath::Power(0.5*multhist->GetBinError(i,j),2)));
486           }
487           else {
488             dNdetadphiHistogramSPD->SetBinContent(i,j,multhist->GetBinContent(i,j));
489             dNdetadphiHistogramSPD->SetBinError(i,j,multhist->GetBinError(i,j));
490           }
491           
492           
493           
494
495         
496         }
497       }
498       delete multhist;
499     }
500   }
501 }
502 //_____________________________________________________________________
503 //
504 //
505 // EOF