Various upgrades and bug fixes etc
[u/mrichter/AliRoot.git] / PWG2 / FORWARD / analysis / AliFMDDndeta.cxx
1
2 #include "AliFMDDndeta.h"
3 #include "TFile.h"
4 #include "AliLog.h"
5 #include "TH1.h"
6 #include "AliFMDAnaParameters.h"
7 #include "AliFMDAnaCalibSharingEfficiency.h"
8 #include "TStyle.h"
9 //#include "TObjArray.h"
10 #include "TCanvas.h"
11 #include "TLine.h"
12 #include "TGraphErrors.h"
13 #include "TGraphAsymmErrors.h"
14 #include "TPad.h"
15 #include "iostream"
16 #include "TH3.h"
17 #include "TMath.h"
18 #include "TProfile.h"
19 #include "TProfile2D.h"
20 #include "TProfile3D.h"
21 #include "TLegend.h"
22 #include "TPad.h"
23 #include "TLatex.h"
24 #include "TStyle.h"
25 #include "TF1.h"
26
27 #define SMALLNUMBER 0.0001
28
29 ClassImp(AliFMDDndeta)
30 //_____________________________________________________________________
31
32 AliFMDDndeta::AliFMDDndeta() 
33 : TObject(),
34   fList(0),
35   fMultList(),
36   fNbinsToCut(2),
37   fVtxCut1(-10),
38   fVtxCut2(10),
39   fIsInit(kFALSE),
40   fIsGenerated(),
41   fPrimEvents(),
42   fEvents(),
43   fPrimdNdeta(),
44   fDrawAll(kFALSE)
45 {
46   //AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
47   /* fDataObject = new TProfile3D("dataObject","dataObject",
48                                pars->GetNetaBins(),-6,6,
49                                20,0,2*TMath::Pi(),
50                                pars->GetNvtxBins(),-0.5,pars->GetNvtxBins()-0.5);
51   fDataObject->SetXTitle("#eta");
52   fDataObject->SetYTitle("#varphi [radians]");
53   fDataObject->SetZTitle("v_{z} [cm]");*/
54   
55   fAnalysisNames[0] = "Hits";
56   fAnalysisNames[1] = "HitsTrVtx";
57   fAnalysisNames[2] = "dNdeta";
58   fAnalysisNames[3] = "dNdetaTrVtx";
59   fAnalysisNames[4] = "dNdetaNSD";
60   
61   for(Int_t i=0; i<5;i++) 
62     fMultList[i] = new TList();
63 }
64 //_____________________________________________________________________
65 void AliFMDDndeta::SetNames(Analysis what) {
66   // Set names of histograms from analysis
67   
68   switch(what) {
69   case kHits :
70     fPrimEvents.Form("nMCEventsNoCuts"); //was nMCEvents
71     fEvents.Form("nEvents");
72     fPrimdNdeta.Form("hMultvsEtaNoCuts");
73     break;
74   case kHitsTrVtx :
75     fPrimEvents.Form("nMCEvents"); 
76     fEvents.Form("nEvents");
77     fPrimdNdeta.Form("hMultvsEtaNoCuts");
78     break;
79   case kMult :
80     fPrimEvents.Form("nMCEventsNoCuts"); 
81     fEvents.Form("nEvents");
82     fPrimdNdeta.Form("hMultvsEtaNoCuts");
83     break;
84   case kMultTrVtx :
85     fPrimEvents.Form("nMCEvents"); 
86     fEvents.Form("nEvents");
87     fPrimdNdeta.Form("hMultvsEta");
88     break;
89   case kMultNSD :
90     fPrimEvents.Form("nMCEventsNSDNoCuts"); 
91     fEvents.Form("nNSDEvents");
92     fPrimdNdeta.Form("hMultvsEtaNSDNoCuts");
93     break;
94     
95   default:
96     break;
97   }
98 }
99 //_____________________________________________________________________
100 const char* AliFMDDndeta::GetAnalysisName(Analysis what, UShort_t det, Char_t ring, Int_t vtxbin) {
101   //Get names of histograms
102   const char* name = "";
103   
104   switch(what) {
105   case kHits :
106     name = Form("hits_NoCuts_FMD%d%c_vtxbin%d_proj",det,ring,vtxbin); //NoCuts added
107     break;
108   case kHitsTrVtx :
109     name = Form("hits_FMD%d%c_vtxbin%d_proj",det,ring,vtxbin); 
110     break;
111   case kMult :
112     name = Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ring,vtxbin);
113     break;
114   case kMultTrVtx :
115     name = Form("dNdeta_FMD%d%c_TrVtx_vtxbin%d_proj",det,ring,vtxbin);
116     break;
117   case kMultNSD :
118     name = Form("dNdetaNSD_FMD%d%c_vtxbin%d_proj",det,ring,vtxbin);
119     break;
120     
121   default:
122     break;
123   }
124   //std::cout<<name.Data()<<std::endl;
125   return name;
126 }
127 //_____________________________________________________________________
128 const char* AliFMDDndeta::GetPrimName(Analysis what, UShort_t det, Char_t ring, Int_t vtxbin) {
129   //Get names of primaries
130   const char* name = "";
131   
132   switch(what) {
133   case kHits :
134     name = Form("hMCHits_nocuts_FMD%d%c_vtxbin%d",det,ring,vtxbin); //nocuts added
135     break;
136   case kHitsTrVtx :
137     name = Form("hMCHits_FMD%d%c_vtxbin%d",det,ring,vtxbin);
138     break;
139   case kMult :
140     name = Form("primmult_vtxbin%d",vtxbin);
141     break;
142   case kMultTrVtx :
143     name = Form("primmult_NoCuts_vtxbin%d",vtxbin);
144     break;
145   case kMultNSD :
146     name = Form("primmult_NSD_vtxbin%d",vtxbin);
147     break;
148   default:
149     break;
150   }
151   return name;
152 }
153 //_____________________________________________________________________
154 void AliFMDDndeta::GenerateHits(Analysis what) {
155   //Generate the hits distributions
156   AliFMDAnaParameters* pars =  AliFMDAnaParameters::Instance();
157   TH2F* hTmp         =  pars->GetBackgroundCorrection(1,'I',0);
158   Int_t nVertexBins  =  pars->GetNvtxBins();
159   Int_t nEtaBins     =  hTmp->GetNbinsX();
160   Float_t etaMin     =  hTmp->GetXaxis()->GetXmin(); 
161   Float_t etaMax     =  hTmp->GetXaxis()->GetXmax();
162
163   for(Int_t i = 0; i<nVertexBins;i++) { // 
164     TH1F* hHits = new TH1F(Form("hMCHits_vtxbin%d_%s",i,fAnalysisNames[what]),Form("hHits_vtxbin%d_%s",i,fAnalysisNames[what]),nEtaBins,etaMin,etaMax);
165     hHits->Sumw2();
166     fMultList[what]->Add(hHits);
167   }
168   
169   for(Int_t det = 1; det<=3; det++) {
170     Int_t maxRing = (det == 1 ? 0 : 1);
171     for(Int_t ring = 0;ring<=maxRing;ring++) {
172       
173       Char_t ringChar = (ring == 0 ? 'I' : 'O');
174       for(Int_t v=0; v< nVertexBins; v++) {
175         TH1F* hits = (TH1F*)fList->FindObject(GetPrimName(what,det,ringChar,v));
176         
177         Int_t nNonZero = 0, nNonZeroInData = 0;
178         
179         //removing edges
180         for(Int_t i =1;i<=hits->GetNbinsX();i++) {
181           if(hits->GetBinContent(i) !=0)
182             nNonZero++;
183         }
184
185         TH1F* sumMultHist = (TH1F*)fMultList[what]->FindObject(Form("hMCHits_vtxbin%d_%s",v,fAnalysisNames[what]));
186         for(Int_t i =1;i<=sumMultHist->GetNbinsX();i++) {
187           if(hits->GetBinContent(i) == 0 ) continue;
188           
189           nNonZeroInData++;
190           
191           if(nNonZeroInData <=fNbinsToCut || nNonZeroInData > (nNonZero - fNbinsToCut)) {
192             continue;
193           }
194           if(sumMultHist->GetBinContent(i) < SMALLNUMBER && TMath::Abs(hits->GetBinContent(i)) > 0){
195               sumMultHist->SetBinContent(i,hits->GetBinContent(i));
196               sumMultHist->SetBinError(i,hits->GetBinError(i));
197             }
198           else {
199               
200               
201             Float_t sumofw = (1/TMath::Power(hits->GetBinError(i),2)) + (1/TMath::Power(sumMultHist->GetBinError(i),2));
202             Float_t wav =  (hits->GetBinContent(i)*(1/TMath::Power(hits->GetBinError(i),2)) + sumMultHist->GetBinContent(i)*(1/TMath::Power(sumMultHist->GetBinError(i),2)))/sumofw;
203             Float_t error = 1/TMath::Sqrt(sumofw);
204             sumMultHist->SetBinContent(i, wav);
205             sumMultHist->SetBinError(i, error);
206           }
207         }
208       }
209     }
210   }
211 }
212
213 //_____________________________________________________________________
214 void AliFMDDndeta::Init(const Char_t* filename) {
215   //Initialize everything from file
216   TFile* fin = TFile::Open(filename);
217   
218   if(!fin) {
219     AliWarning("No file - exiting !");
220     return;
221   }
222   AliFMDAnaParameters* pars =  AliFMDAnaParameters::Instance();
223   //pars->Init();
224   
225   TList* list = (TList*)fin->Get(Form("%s/BackgroundCorrected",pars->GetDndetaAnalysisName()));
226   
227   if(!list) //an old file ? Perhaps...
228     list = (TList*)fin->Get("BackgroundCorrected");
229   
230   Init(list);
231   
232 }
233 //_____________________________________________________________________
234 void AliFMDDndeta::Init(TList* list) {
235   //Initialize everything from TList
236     
237   if(!list) {
238     AliWarning("No list - exiting !");
239     return;
240   }
241   
242   fList = (TList*)list->Clone("inputList");
243     
244   fIsGenerated[kHits]      = kFALSE;
245   fIsGenerated[kMult]      = kFALSE;  
246   fIsGenerated[kMultTrVtx] = kFALSE;
247   fIsGenerated[kHitsTrVtx] = kFALSE;
248   fIsGenerated[kMultNSD]   = kFALSE;
249   
250   fIsInit = kTRUE;
251 }
252 //_____________________________________________________________________
253 void AliFMDDndeta::GenerateMult(Analysis what) {
254   //Generate dNdeta
255   
256   if(!fIsInit) {
257     AliWarning("Not initialised - call Init to remedy");
258     return;
259   }
260   
261   if(fIsGenerated[what]) {
262     AliWarning("Already generated - have a look at the results!");
263     return;
264   }
265   else fIsGenerated[what] = kTRUE;
266   
267   SetNames(what);
268   
269   if(what == kHits || what == kHitsTrVtx)
270     GenerateHits(what);
271   
272   TH1I* hMCEvents = (TH1I*)fList->FindObject(fPrimEvents.Data());
273   
274   AliFMDAnaParameters* pars =  AliFMDAnaParameters::Instance();
275   
276   TH2F* hTmp         =  pars->GetBackgroundCorrection(1,'I',0);
277   Int_t nVertexBins  =  pars->GetNvtxBins();
278   Int_t nEtaBins     =  hTmp->GetNbinsX();
279   Float_t etaMin     =  hTmp->GetXaxis()->GetXmin(); 
280   Float_t etaMax     =  hTmp->GetXaxis()->GetXmax();
281   
282   for(Int_t i = 0; i<nVertexBins;i++) {
283     TH1F* hMult = new TH1F(Form("hMult_vtxbin%d_%s",i,fAnalysisNames[what]),Form("hMult_vtxbin%d_%s",i,fAnalysisNames[what]),nEtaBins,etaMin,etaMax);
284     hMult->Sumw2();
285     fMultList[what]->Add(hMult);
286   }
287   
288   for(Int_t det = 1; det<=3;det++) {
289     Int_t maxRing = (det == 1 ? 0 : 1);
290     for(Int_t iring = 0; iring<=maxRing; iring++) {
291       Char_t ringChar = (iring == 0 ? 'I' : 'O');
292       //Int_t nsec  = (iring == 0 ? 20 : 40);
293       TH1F* hRingMult= new TH1F(Form("hRingMult_FMD%d%c_%s",det,ringChar,fAnalysisNames[what]),Form("hRingMult_FMD%d%c_%s",det,ringChar,fAnalysisNames[what]),nEtaBins,etaMin,etaMax);
294       fMultList[what]->Add(hRingMult);
295       //      TProfile* phiprofile     = new TProfile(Form("dNdphiFMD%d%c",det,ringChar), Form("dNdphiFMD%d%c;#Phi",det,ringChar), nsec , 0, 2*TMath::Pi());
296       //fMultList[what]->Add(phiprofile);
297     }
298   }
299   TH1I*  hEvents  = (TH1I*)fList->FindObject(fEvents.Data());
300   TH1F*  hPrimVtx = 0;
301   
302         
303     
304   for(Int_t det = 1; det<=3; det++) {
305     Int_t maxRing = (det == 1 ? 0 : 1);
306     for(Int_t ring = 0;ring<=maxRing;ring++) {
307       
308       Char_t ringChar = (ring == 0 ? 'I' : 'O');
309       for(Int_t v=0; v< nVertexBins; v++) {
310         if(det == 1) {
311           if(what == kHits || what == kHitsTrVtx)
312             hPrimVtx = (TH1F*)fMultList[what]->FindObject(Form("hMCHits_vtxbin%d_%s",v,fAnalysisNames[what]));
313           else
314             hPrimVtx = (TH1F*)fList->FindObject(GetPrimName(what,det,ringChar,v));
315           
316
317           Float_t nPrimevents = hMCEvents->GetBinContent(v+1);
318           Float_t xb = hPrimVtx->GetNbinsX();
319           Float_t xr = hPrimVtx->GetXaxis()->GetXmax() - hPrimVtx->GetXaxis()->GetXmin(); 
320           hPrimVtx->Scale(xb / xr );
321           if(nPrimevents > 0)
322             hPrimVtx->Scale(1/nPrimevents);
323           
324         }
325         
326         Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
327         Float_t vtxZ1   = (delta*v) - pars->GetVtxCutZ();
328         Float_t vtxZ2   = (delta*(v+1)) - pars->GetVtxCutZ();
329         
330         if(vtxZ1< fVtxCut1 || vtxZ2 >fVtxCut2)
331           continue;
332         Float_t nEvents = (Float_t)hEvents->GetBinContent(v+1);
333         
334         //TH1F* multhistproj = (TH1F*)fList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d_proj",det,ringChar,v));
335         
336
337         
338         
339         
340         
341         TH1F* multhistproj = (TH1F*)fList->FindObject(GetAnalysisName(what,det,ringChar,v));
342         if(nEvents)
343           multhistproj->Scale(1/nEvents);
344         Float_t xnBins = multhistproj->GetNbinsX();
345         Float_t xrange = multhistproj->GetXaxis()->GetXmax() - multhistproj->GetXaxis()->GetXmin(); 
346         multhistproj->Scale(xnBins / xrange);
347         
348         //TH2F* multhist = (TH2F*)fList->FindObject(Form("dNdeta_FMD%d%c_vtxbin%d",det,ringChar,v));
349
350         
351         //if(nEvents)
352         //  multhist->Scale(1/nEvents);
353         /*      for(Int_t i=1;i<=multhist->GetNbinsX();i++) {
354           for(Int_t j=1;j<=multhist->GetNbinsY();j++) {
355             
356             if (multhist->GetBinContent(i,j) <= 0.0001) continue;
357             //std::cout<<multhist->GetBinContent(i,j)<<std::endl;
358             fDataObject->Fill(multhist->GetXaxis()->GetBinCenter(i),
359                               multhist->GetYaxis()->GetBinCenter(j),
360                               v,
361                               multhist->GetBinContent(i,j), multhist->GetBinError(i,j));
362             //fDataObject->SetBinContent(i,j,v,multhist->GetBinContent(i,j));
363             //fDataObject->SetBinError(i,j,v,multhist->GetBinError(i,j));
364           } 
365         }
366         */      
367         //if(nEvents)
368         //  multhist->Scale(1/nEvents);
369         //if(ringChar == 'O')
370         //  multhist->RebinY(2);
371         //removing edges
372         
373         Int_t nNonZero = 0, nNonZeroInData = 0;
374         
375         for(Int_t i =1;i<=multhistproj->GetNbinsX();i++) {
376           if(multhistproj->GetBinContent(i) !=0)
377             nNonZero++;
378         }
379         Int_t nBinsOld = fNbinsToCut;
380         //if(det == 1 && ringChar =='I') {
381         //  fNbinsToCut = 0;
382         //      }
383         TH1F* hRingMult = (TH1F*)fMultList[what]->FindObject(Form("hRingMult_FMD%d%c_%s",det,ringChar,fAnalysisNames[what]));
384         
385         for(Int_t i=1; i<=hRingMult->GetNbinsX(); i++) {
386           if(multhistproj->GetBinContent(i)!=0) {
387             nNonZeroInData++;
388             if(nNonZeroInData <=fNbinsToCut || nNonZeroInData > (nNonZero - fNbinsToCut)) {
389               continue;
390             }
391             Float_t oldweight = 0;
392             Float_t oldwav    = 0;
393             if(hRingMult->GetBinError(i)>0) {
394               oldweight = 1/TMath::Power(hRingMult->GetBinError(i),2);
395               oldwav    = oldweight*hRingMult->GetBinContent(i);
396             }
397             Float_t  weight = 1/TMath::Power(multhistproj->GetBinError(i),2);
398             Float_t  wav    = oldwav + weight*multhistproj->GetBinContent(i);
399             Float_t  sumofw = oldweight + weight;
400             if(sumofw) {
401               Float_t error = 1/TMath::Sqrt(sumofw);
402               
403               hRingMult->SetBinContent(i,wav/sumofw);
404               hRingMult->SetBinError(i,error);
405               
406             }
407           }
408         }
409         nNonZeroInData = 0;
410         
411         TH1F* sumMultHist = (TH1F*)fMultList[what]->FindObject(Form("hMult_vtxbin%d_%s",v,fAnalysisNames[what]));
412         
413         for(Int_t i =1;i<=sumMultHist->GetNbinsX();i++) {
414           if(multhistproj->GetBinContent(i) != 0 ) {      
415             nNonZeroInData++;
416           
417             if(nNonZeroInData <=fNbinsToCut || nNonZeroInData > (nNonZero - fNbinsToCut)) {
418               continue;
419             }
420             if(sumMultHist->GetBinContent(i) < SMALLNUMBER && TMath::Abs(multhistproj->GetBinContent(i)) > 0){
421               sumMultHist->SetBinContent(i,multhistproj->GetBinContent(i));
422               sumMultHist->SetBinError(i,multhistproj->GetBinError(i));
423             }
424             else {
425               
426               
427               Float_t sumofw = (1/TMath::Power(multhistproj->GetBinError(i),2)) + (1/TMath::Power(sumMultHist->GetBinError(i),2));
428               Float_t wav =  (multhistproj->GetBinContent(i)*(1/TMath::Power(multhistproj->GetBinError(i),2)) + sumMultHist->GetBinContent(i)*(1/TMath::Power(sumMultHist->GetBinError(i),2)))/sumofw;
429               Float_t error = 1/TMath::Sqrt(sumofw);
430               sumMultHist->SetBinContent(i, wav);
431               sumMultHist->SetBinError(i, error);
432             }
433             
434         
435           }
436         }
437         fNbinsToCut = nBinsOld;
438       }
439       
440     }
441   }
442   
443   TH1F* sumMult  = new TH1F(Form("dNdeta_%s",fAnalysisNames[what]),Form("dNdeta_%s",fAnalysisNames[what]),nEtaBins,etaMin,etaMax);
444   sumMult->Sumw2();
445   fMultList[what]->Add(sumMult);
446   TH1F* primMult  = new TH1F(Form("primary_%s",fAnalysisNames[what]),Form("primary_%s",fAnalysisNames[what]),nEtaBins,etaMin,etaMax);
447   primMult->Sumw2();
448   fMultList[what]->Add(primMult);
449
450   Float_t wav  = 0, sumofw=0, weight = 0, primwav  = 0, primsumofw=0, primweight = 0;//, etaofbin = 0;
451   for(Int_t i =1; i<=sumMult->GetNbinsX();i++) {
452     
453     for(Int_t v = 0; v<nVertexBins;v++) {
454       if(what == kHits || what == kHitsTrVtx)
455         hPrimVtx = (TH1F*)fMultList[what]->FindObject(Form("hMCHits_vtxbin%d_%s",v,fAnalysisNames[what]));
456       else
457         hPrimVtx = (TH1F*)fList->FindObject(GetPrimName(what,0,0,v));
458       TH1F* sumMultHist = (TH1F*)fMultList[what]->FindObject(Form("hMult_vtxbin%d_%s",v,fAnalysisNames[what]));
459       //etaofbin += sumMultHist->GetBinCenter(i);
460       // if( hPrimVtx->GetBinContent(i)!=0 && hPrimVtx->GetBinError(i)>0.001 && sumMultHist->GetBinContent(i)!=0) {
461       if( TMath::Abs(hPrimVtx->GetBinContent(i)) > 0) {
462         primweight = 1/TMath::Power(hPrimVtx->GetBinError(i),2);
463         primwav    = primwav + primweight*hPrimVtx->GetBinContent(i);
464         primsumofw = primsumofw + primweight;
465       }
466       if(TMath::Abs(sumMultHist->GetBinContent(i)) > 0) {
467         
468         weight = 1/TMath::Power(sumMultHist->GetBinError(i),2);
469         wav    = wav + weight*sumMultHist->GetBinContent(i);
470         sumofw = sumofw + weight;
471       }
472       
473     }
474     if( primsumofw !=0 ) {// && sumofw !=0) {
475       Float_t primerror = 1/TMath::Sqrt(primsumofw);
476       
477       primMult->SetBinContent(i,primwav/primsumofw);
478       primMult->SetBinError(i,primerror);
479     }
480     else {
481       primMult->SetBinContent(i,0);
482       primMult->SetBinError(i,0);
483     }
484     
485     if(sumofw) {
486       
487       Float_t error = 1/TMath::Sqrt(sumofw);
488       
489       sumMult->SetBinContent(i,wav/sumofw);
490       sumMult->SetBinError(i,error);
491     }
492     
493     
494     wav = 0;
495     sumofw = 0;
496     weight = 0;
497     primwav = 0;
498     primsumofw = 0;
499     primweight = 0;
500     
501   }
502   
503 }
504 //_____________________________________________________________________
505 void AliFMDDndeta::DrawDndeta(Analysis what, Int_t rebin, Bool_t realdata, TString pythiafile) {
506   //Draw everything
507   if(!fIsInit) {
508     AliWarning("Not initialised - call Init to remedy");
509     return;
510   }
511
512   if(!fIsGenerated[what]) {
513     AliWarning("Not generated - generate first then draw!");
514     return;
515   }
516   AliFMDAnaParameters* pars =  AliFMDAnaParameters::Instance();
517   SetNames(what);
518   //UA5 data NSD
519   Float_t x[19] = {0.125,0.375,0.625,0.875,1.125,1.375,1.625,1.875,2.125,2.375,
520                    2.625,2.875,3.125,3.375,3.625,3.875,4.125,4.375,4.625};
521   Float_t x2[19] = {-0.125,-0.375,-0.625,-0.875,-1.125,-1.375,-1.625,-1.875,
522                     -2.125,-2.375,-2.625,-2.875,-3.125,-3.375,-3.625,-3.875,
523                     -4.125,-4.375,-4.625};
524   
525   
526   Float_t y[19] = {3.48, 3.38, 3.52, 3.68, 3.71, 3.86, 3.76, 3.66, 3.72 ,3.69,
527                    3.56, 3.41, 3.15, 3.09, 2.74, 2.73, 2.32, 1.99, 1.69};
528   
529   Float_t ey[19] = {0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,0.07,
530                     0.07,0.07,0.08,0.08,0.09,0.09,0.1,0.13};
531     
532   TGraphErrors* graph = new TGraphErrors(19,x,y,NULL,ey);
533   TGraphErrors* graph2 = new TGraphErrors(19,x2,y,NULL,ey);
534   //UA5 data INEL
535   Float_t xinel[19] = {0.125, 0.375, 0.625, 0.875, 1.125, 1.375, 1.625, 1.875, 2.125, 2.375, 2.625, 2.875, 3.125, 3.375, 3.625, 3.875, 4.125, 4.375, 4.625};
536   Float_t xinel2[19] = {-0.125, -0.375, -0.625, -0.875, -1.125, -1.375, -1.625, -1.875, -2.125, -2.375, -2.625, -2.875, -3.125, -3.375, -3.625, -3.875, -4.125, -4.375, -4.625};
537   Float_t yinel[19] = {3.14, 3.04, 3.17, 3.33, 3.33, 3.53, 3.46, 3.41, 3.45, 3.39, 3.07, 3.07, 2.93, 2.93, 2.55, 2.48, 2.18, 1.91, 1.52};
538   Float_t eyinel[19] = {0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.07, 0.08, 0.08, 0.09, 0.09, 0.1, 0.13};
539  TGraphErrors* graphinel = new TGraphErrors(19,xinel,yinel,NULL,eyinel);
540  TGraphErrors* graphinel2 = new TGraphErrors(19,xinel2,yinel,NULL,eyinel);
541   
542   graph->SetMarkerStyle(22);
543   graph->SetMarkerColor(kBlue);
544   graphinel->SetMarkerStyle(20);
545   graphinel->SetMarkerColor(kGreen);
546   graphinel2->SetMarkerStyle(24);
547   graphinel2->SetMarkerColor(kGreen);
548   //graph->Draw("P");
549   graph2->SetMarkerStyle(26);
550   graph2->SetMarkerColor(kBlue);
551   graphinel->GetHistogram()->SetStats(kFALSE);
552   graphinel2->GetHistogram()->SetStats(kFALSE);
553   
554   
555   //End UA5 data
556   
557   //Published ALICE data
558   // Plot: p7742_d1x1y1
559   
560   // INEL
561   
562   TGraphAsymmErrors* p7742_d1x1y1 = 0;
563   double p7742_d1x1y1_xval[] = { -1.3, -1.1, -0.9, -0.7, -0.5, -0.3, -0.1, 0.1, 0.3, 
564     0.5, 0.7, 0.9, 1.1, 1.3 };
565   double p7742_d1x1y1_xerrminus[] = { 0.09999999999999987, 0.09999999999999987, 0.09999999999999998, 0.10000000000000009, 0.09999999999999998, 0.10000000000000003, 0.1, 0.1, 0.09999999999999998, 
566     0.09999999999999998, 0.09999999999999998, 0.09999999999999998, 0.10000000000000009, 0.10000000000000009 };
567   double p7742_d1x1y1_xerrplus[] = { 0.10000000000000009, 0.10000000000000009, 0.09999999999999998, 0.09999999999999998, 0.09999999999999998, 0.09999999999999998, 0.1, 0.1, 0.10000000000000003, 
568     0.09999999999999998, 0.10000000000000009, 0.09999999999999998, 0.09999999999999987, 0.09999999999999987 };
569   double p7742_d1x1y1_yval[] = { 3.28, 3.28, 3.22, 3.12, 3.06, 3.02, 2.98, 3.02, 3.02, 
570     3.05, 3.15, 3.21, 3.26, 3.33 };
571   double p7742_d1x1y1_yerrminus[] = { 0.06324555320336758, 0.06324555320336758, 0.06324555320336758, 0.06324555320336758, 0.06324555320336758, 0.05385164807134505, 0.05385164807134505, 0.05385164807134505, 0.05385164807134505, 
572     0.06324555320336758, 0.06324555320336758, 0.06324555320336758, 0.06324555320336758, 0.06324555320336758 };
573   double p7742_d1x1y1_yerrplus[] = { 0.08246211251235322, 0.08246211251235322, 0.08246211251235322, 0.08246211251235322, 0.08246211251235322, 0.08246211251235322, 0.07280109889280519, 0.08246211251235322, 0.08246211251235322, 
574     0.08246211251235322, 0.08246211251235322, 0.08246211251235322, 0.08246211251235322, 0.08246211251235322 };
575   int p7742_d1x1y1_numpoints = 14;
576   p7742_d1x1y1 = new TGraphAsymmErrors(p7742_d1x1y1_numpoints, p7742_d1x1y1_xval, p7742_d1x1y1_yval, p7742_d1x1y1_xerrminus, p7742_d1x1y1_xerrplus, p7742_d1x1y1_yerrminus, p7742_d1x1y1_yerrplus);
577   p7742_d1x1y1->SetName("/HepData/7742/d1x1y1");
578   p7742_d1x1y1->SetTitle("/HepData/7742/d1x1y1");
579   p7742_d1x1y1->SetMarkerStyle(21);
580   p7742_d1x1y1->SetMarkerColor(kRed);
581   // p7742_d1x1y1->Draw("Psame");
582   
583   
584   //NSD
585   TGraphAsymmErrors*  p7742_d2x1y1 = 0;
586  double p7742_d2x1y1_xval[] = { -1.3, -1.1, -0.9, -0.7, -0.5, -0.3, -0.1, 0.1, 0.3, 
587     0.5, 0.7, 0.9, 1.1, 1.3 };
588   double p7742_d2x1y1_xerrminus[] = { 0.09999999999999987, 0.09999999999999987, 0.09999999999999998, 0.10000000000000009, 0.09999999999999998, 0.10000000000000003, 0.1, 0.1, 0.09999999999999998, 
589     0.09999999999999998, 0.09999999999999998, 0.09999999999999998, 0.10000000000000009, 0.10000000000000009 };
590   double p7742_d2x1y1_xerrplus[] = { 0.10000000000000009, 0.10000000000000009, 0.09999999999999998, 0.09999999999999998, 0.09999999999999998, 0.09999999999999998, 0.1, 0.1, 0.10000000000000003, 
591     0.09999999999999998, 0.10000000000000009, 0.09999999999999998, 0.09999999999999987, 0.09999999999999987 };
592   double p7742_d2x1y1_yval[] = { 3.9, 3.89, 3.81, 3.7, 3.64, 3.59, 3.53, 3.58, 3.59, 
593     3.61, 3.74, 3.8, 3.87, 3.95 };
594   double p7742_d2x1y1_yerrminus[] = { 0.13341664064126335, 0.13152946437965907, 0.13152946437965907, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644, 
595     0.1216552506059644, 0.1216552506059644, 0.13152946437965907, 0.13152946437965907, 0.13341664064126335 };
596   double p7742_d2x1y1_yerrplus[] = { 0.13341664064126335, 0.13152946437965907, 0.13152946437965907, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644, 0.1216552506059644, 
597     0.1216552506059644, 0.1216552506059644, 0.13152946437965907, 0.13152946437965907, 0.13341664064126335 };
598   int p7742_d2x1y1_numpoints = 14;
599
600   p7742_d2x1y1 = new TGraphAsymmErrors(p7742_d2x1y1_numpoints, p7742_d2x1y1_xval, p7742_d2x1y1_yval, p7742_d2x1y1_xerrminus, p7742_d2x1y1_xerrplus, p7742_d2x1y1_yerrminus, p7742_d2x1y1_yerrplus);
601   p7742_d2x1y1->SetName("/HepData/7742/d2x1y1");
602   p7742_d2x1y1->SetTitle("/HepData/7742/d2x1y1");
603   p7742_d2x1y1->SetMarkerStyle(21);
604   p7742_d2x1y1->SetMarkerColor(kRed);
605   //p7742_d2x1y1.Draw("AP");
606
607
608
609
610
611
612   //End official data
613   
614   // CMS published NSD data
615   
616   TGraphAsymmErrors* p7743_d8x1y1 = 0;
617   double p7743_d8x1y1_xval[] = { -2.25, -1.75, -1.25, -0.75, -0.25, 0.25, 0.75, 1.25, 1.75, 
618     2.25 };
619   double p7743_d8x1y1_xerrminus[] = { 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25 };
620   double p7743_d8x1y1_xerrplus[] = { 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25, 0.25 };
621   double p7743_d8x1y1_yval[] = { 3.6, 3.73, 3.62, 3.54, 3.48, 3.48, 3.54, 3.62, 3.73,  3.6 };
622   double p7743_d8x1y1_yerrminus[] = { 0.13, 0.14, 0.13, 0.13, 0.13, 0.13, 0.13, 0.13, 0.14,0.13 };
623   double p7743_d8x1y1_yerrplus[] = { 0.13, 0.14, 0.13, 0.13, 0.13, 0.13, 0.13, 0.13, 0.14, 0.13 };
624   int p7743_d8x1y1_numpoints = 10;
625   p7743_d8x1y1 = new TGraphAsymmErrors(p7743_d8x1y1_numpoints, p7743_d8x1y1_xval, p7743_d8x1y1_yval, p7743_d8x1y1_xerrminus, p7743_d8x1y1_xerrplus, p7743_d8x1y1_yerrminus, p7743_d8x1y1_yerrplus);
626
627   p7743_d8x1y1->SetMarkerStyle(20);
628   p7743_d8x1y1->SetMarkerColor(kBlack);
629   
630   //End CMS
631   
632   
633   TH1I* hEvents   = (TH1I*)fList->FindObject(fEvents.Data());
634   
635   TH1I* hMCEvents = (TH1I*)fList->FindObject(fPrimEvents.Data());
636   
637    //SPD part
638   TH1D* hSPDana = (TH1D*)fList->FindObject(Form("dNdeta_SPD_vtxbin%d_proj",5));
639   
640   for(Int_t nn = 0; nn < pars->GetNvtxBins() ; nn++) {
641     TH1D* hSPDanalysis = (TH1D*)fList->FindObject(Form("dNdeta_SPD_vtxbin%d_proj",nn));
642     TH1D* hSPDanalysisTrVtx = (TH1D*)fList->FindObject(Form("dNdetaTrVtx_SPD_vtxbin%d_proj",nn));
643     TH1D* hSPDanalysisNSD = (TH1D*)fList->FindObject(Form("dNdetaNSD_SPD_vtxbin%d_proj",nn));
644     
645     Float_t nEventsSPD = (Float_t)hEvents->GetBinContent(nn+1);
646     if(!nEventsSPD) continue; 
647     hSPDanalysis->Scale(1/nEventsSPD);
648     hSPDanalysisTrVtx->Scale(1/nEventsSPD);
649     hSPDanalysisNSD->Scale(1/nEventsSPD);
650   }
651   
652   TH1F* hSPDcombi = new TH1F("SPDcombi","SPDcombi",hSPDana->GetNbinsX(),hSPDana->GetXaxis()->GetXmin(),hSPDana->GetXaxis()->GetXmax());
653   TH1D* hSPDanalysis = 0;
654   for(Int_t kk = 1; kk <=hSPDana->GetNbinsX(); kk++) {
655     Float_t weight = 0, wav=0,sumofw = 0;
656     for(Int_t nn = 0; nn < pars->GetNvtxBins() ; nn++) {
657       if(what == kMult)
658         hSPDanalysis = (TH1D*)fList->FindObject(Form("dNdeta_SPD_vtxbin%d_proj",nn));
659       else if(what == kMultNSD)
660         hSPDanalysis = (TH1D*)fList->FindObject(Form("dNdetaNSD_SPD_vtxbin%d_proj",nn));
661       else if(what == kMultTrVtx)
662         hSPDanalysis = (TH1D*)fList->FindObject(Form("dNdetaTrVtx_SPD_vtxbin%d_proj",nn));
663       else continue;
664       
665       if(TMath::Abs(hSPDanalysis->GetBinCenter(kk)) > 2)
666         continue;
667       Float_t mult = hSPDanalysis->GetBinContent(kk);
668       Float_t error = hSPDanalysis->GetBinError(kk);
669       
670       if(mult > 0 && hSPDanalysis->GetBinContent(kk-1) < SMALLNUMBER)
671         continue;
672       if(mult > 0 && hSPDanalysis->GetBinContent(kk-2) < SMALLNUMBER)
673         continue;
674       
675       if(mult > 0 && hSPDanalysis->GetBinContent(kk+1) < SMALLNUMBER)
676         continue;
677       if(mult > 0 && hSPDanalysis->GetBinContent(kk+2) < SMALLNUMBER)
678         continue;
679       
680       if(mult > 0) {
681         weight = 1/TMath::Power(error,2);
682         wav    = wav + weight*mult;
683         sumofw = sumofw + weight;
684       }
685       
686     }
687     
688     if(sumofw && wav) {
689       Float_t errorTotal = 1/TMath::Sqrt(sumofw);
690       hSPDcombi->SetBinContent(kk,wav/sumofw);
691       hSPDcombi->SetBinError(kk,errorTotal);
692     }
693   }
694     
695   
696   Float_t xb1 = hSPDcombi->GetNbinsX();
697   Float_t xr1 = hSPDcombi->GetXaxis()->GetXmax() - hSPDcombi->GetXaxis()->GetXmin(); 
698   hSPDcombi->Scale(xb1 / xr1);
699   //hSPDcombi->Rebin(rebin);
700   //hSPDcombi->Scale(1/(Float_t)rebin);
701   //RebinHistogram(hSPDcombi,rebin);
702   hSPDcombi->SetMarkerStyle(29);
703   hSPDcombi->SetMarkerColor(kBlue);
704   //hSPDcombi->DrawCopy("same");
705   
706   
707   TCanvas* c1 = new TCanvas("cvertexbins","Overlaps",1400,800);
708   c1->Divide(5,2);
709   TCanvas* cCorrections = 0;
710   TCanvas* cCorrectionsPhi = 0;
711   if(fDrawAll) {
712     cCorrections = new TCanvas("corrections","Correction vs Eta",1400,800);
713     cCorrections->Divide(5,2);
714     
715     cCorrectionsPhi = new TCanvas("correctionsPhi","Correction vs Phi",1400,800);
716     cCorrectionsPhi->Divide(5,2);
717   }
718   //TCanvas* cphi = new TCanvas("dNdphi","dNdphi",1280,1024);
719   // cphi->Divide(3,2);
720   
721   Int_t nVertexBins  =  pars->GetNvtxBins();
722   
723   
724   //Int_t npadphi = 0;
725   TH1F*  hPrimVtx = 0;
726   for(Int_t det = 1; det<=3; det++) {
727     Int_t maxRing = (det == 1 ? 0 : 1);
728     for(Int_t ring = 0;ring<=maxRing;ring++) {
729       //  npadphi++;
730       
731       for(Int_t v=0; v< nVertexBins; v++) {
732         Char_t ringChar = (ring == 0 ? 'I' : 'O');
733         Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
734         Float_t vtxZ1   = (delta*v) - pars->GetVtxCutZ();
735         Float_t vtxZ2   = (delta*(v+1)) - pars->GetVtxCutZ();
736           
737         if(vtxZ1<fVtxCut1 || vtxZ2 >fVtxCut2)
738           continue;
739         Float_t nEvents = (Float_t)hEvents->GetBinContent(v+1);
740         if(fDrawAll) {
741         TH2F* hBgCor  = pars->GetBackgroundCorrection(det,ringChar,v);
742         if(hBgCor)  {
743           TH1D* hBgProj = hBgCor->ProjectionX(Form("hBgProj_FMD%d%c_vtxbin%d",det,ringChar,v));
744           TH1D* hBgProjPhi = hBgCor->ProjectionY(Form("hBgProjPhi_FMD%d%c_vtxbin%d",det,ringChar,v));
745           
746           hBgProjPhi->SetName(Form("FMD%d%c",det,ringChar));
747           hBgProjPhi->SetTitle("");//Form("FMD%d",det));
748           hBgProj->SetTitle("");
749           //Float_t scalefactor = (hBgProj->GetXaxis()->GetXmax() - hBgProj->GetXaxis()->GetXmin()) / (Float_t)hBgProj->GetNbinsX();
750           Float_t scalefactor = 1/(Float_t)hBgCor->GetNbinsY();
751           
752           Float_t nFilledEtaBins = 0;
753           
754           for(Int_t jj = 1; jj<=hBgProj->GetNbinsX(); jj++) {
755             if(hBgProj->GetBinContent(jj) > 0.01)
756               nFilledEtaBins++;
757           }
758           
759           Float_t scalefactorPhi = 1/nFilledEtaBins;
760           
761           
762           hBgProj->Scale(scalefactor);
763           cCorrections->cd(v+1);
764           
765           hBgProj->SetMarkerColor(det + 2*ring);
766           hBgProj->SetMarkerStyle(19 + det + 2*ring );
767           
768           if((det + 2*ring) == 5)  hBgProj->SetMarkerColor(det + 2*ring+1);
769           if((19 + det + 2*ring ) == 24)  hBgProj->SetMarkerStyle(29 );
770           
771           if(det == 1) {
772             hBgProj->GetYaxis()->SetRangeUser(0,4);
773             hBgProj->SetStats(kFALSE);
774             hBgProj->SetXTitle("#eta");
775             hBgProj->DrawCopy("PE");
776             TLatex* l = new TLatex(0.14,0.92,Form("Vtx range [%.1f, %.1f]",vtxZ1,vtxZ2));
777             l->SetNDC(kTRUE);
778             l->Draw();
779           }
780           else
781             hBgProj->DrawCopy("PEsame");
782         
783           hBgProjPhi->Scale(scalefactorPhi);
784           cCorrectionsPhi->cd(v+1);
785           hBgProjPhi->SetMarkerColor(det + 2*ring);
786           hBgProjPhi->SetMarkerStyle(19 + det + 2*ring );
787           if((det + 2*ring) == 5)  hBgProjPhi->SetMarkerColor(det + 2*ring+1);
788           if((19 + det + 2*ring ) == 24)  hBgProjPhi->SetMarkerStyle(29 );
789           if(det == 1) {
790             hBgProjPhi->GetYaxis()->SetRangeUser(1,5);
791             hBgProjPhi->SetStats(kFALSE);
792             
793             hBgProjPhi->SetXTitle("#Phi");
794             hBgProjPhi->DrawCopy("PE");
795            
796             
797           }
798           else
799             hBgProjPhi->DrawCopy("PEsame");
800           
801         }
802         
803         }
804         
805         TH1F* multhistproj = (TH1F*)fList->FindObject(GetAnalysisName(what,det,ringChar,v));
806                 
807         c1->cd(v+1);
808         
809         if(det==1)  {
810           multhistproj->SetMarkerStyle(20); 
811           multhistproj->SetMarkerColor(1); 
812         }
813         if(det==2 && ringChar=='I') {
814           multhistproj->SetMarkerStyle(21);
815           multhistproj->SetMarkerColor(2); 
816         }
817         if(det==2 && ringChar=='O') {
818           multhistproj->SetMarkerStyle(3);
819           multhistproj->SetMarkerColor(3); 
820         }
821         if(det==3 && ringChar=='I') {
822           multhistproj->SetMarkerStyle(22);
823           multhistproj->SetMarkerColor(4); 
824         }
825         if(det==3 && ringChar=='O') {
826           multhistproj->SetMarkerStyle(23);
827           multhistproj->SetMarkerColor(6); 
828         }
829         
830         
831         
832         if(det == 1) {
833           if(what == kHits || what == kHitsTrVtx)
834             hPrimVtx = (TH1F*)fMultList[what]->FindObject(Form("hMCHits_vtxbin%d_%s",v,fAnalysisNames[what]));
835           else
836             hPrimVtx = (TH1F*)fList->FindObject(GetPrimName(what,det,ringChar,v));
837           // std::cout<<hPrimVtx<<"   "<<kHits<<"   "<<kHitsTrVtx<<"   "<<what<<std::endl;
838           //std::cout<<Form("hMCHits_vtxbin%d_%s",v,fAnalysisNames[what])<<std::endl;
839           hPrimVtx->SetTitle("");
840           TLatex* l = new TLatex(0.14,0.92,Form("Vtx range [%.1f, %.1f], %d events",vtxZ1,vtxZ2,(Int_t)nEvents));
841           l->SetNDC(kTRUE);
842           hPrimVtx->SetFillColor(kGray);
843           hPrimVtx->SetLabelFont(132,"xyz");
844           hPrimVtx->SetStats(0000);
845           hPrimVtx->GetXaxis()->SetTitle("#eta");
846           hPrimVtx->GetYaxis()->SetRangeUser(0,6);
847           hPrimVtx->DrawCopy("E3");
848           l->Draw();
849           multhistproj->DrawCopy("same");
850           
851           TH1D* hSPDanavtxbin = (TH1D*)fList->FindObject(Form("dNdeta_SPD_vtxbin%d_proj",v));
852           hSPDanavtxbin->SetMarkerColor(kBlue);
853           hSPDanavtxbin->SetMarkerStyle(kStar);
854           hSPDanavtxbin->Scale(xb1 / xr1);
855           hSPDanavtxbin->DrawCopy("same");
856           
857           if(what != kMultNSD) {
858             graphinel->Draw("sameP");
859             graphinel2->Draw("sameP");
860           }
861           else
862            {
863              graph->Draw("sameP");
864              graph2->Draw("sameP");
865            }
866         }
867         else
868           multhistproj->DrawCopy("same");
869         
870         
871       }
872     }
873   }
874   
875   //Legends for corrections
876   if(fDrawAll) {
877   for(Int_t v=0; v< nVertexBins; v++) {
878     TPad* pad= (TPad*)cCorrectionsPhi->cd(v+1);
879     pad->BuildLegend(0.15,0.45,0.45,0.9);
880     Double_t delta = 2*pars->GetVtxCutZ()/pars->GetNvtxBins();
881     Float_t vtxZ1   = (delta*v) - pars->GetVtxCutZ();
882     Float_t vtxZ2   = (delta*(v+1)) - pars->GetVtxCutZ();
883     TLatex* l = new TLatex(0.14,0.92,Form("Vtx range [%.1f, %.1f]",vtxZ1,vtxZ2));
884     l->SetNDC(kTRUE);
885     l->Draw();
886   }
887   }
888   for(Int_t v=0; v< nVertexBins; v++) {
889     TH1F* sumMultHist = (TH1F*)fMultList[what]->FindObject(Form("hMult_vtxbin%d_%s",v,fAnalysisNames[what]));
890     c1->cd(v+1);
891     sumMultHist->SetMarkerStyle(25);
892     sumMultHist->DrawCopy("same");
893     
894   }
895   TH1F* primMult = (TH1F*)fMultList[what]->FindObject(Form("primary_%s",fAnalysisNames[what]));
896   TH1F* sumMult  = (TH1F*)fMultList[what]->FindObject(Form("dNdeta_%s",fAnalysisNames[what]));
897   sumMult->SetMarkerStyle(23);
898   sumMult->SetMarkerColor(kRed);
899   for(Int_t det = 1; det<=3;det++) {
900     Int_t maxRing = (det == 1 ? 0 : 1);
901     for(Int_t iring = 0; iring<=maxRing; iring++) {
902       Char_t ringChar = (iring == 0 ? 'I' : 'O');
903       TH1F* hRingMult= (TH1F*)fMultList[what]->FindObject(Form("hRingMult_FMD%d%c_%s",det,ringChar,fAnalysisNames[what]));
904       hRingMult->Add(primMult,-1);
905       hRingMult->Divide(primMult);
906     }
907   }
908   
909  
910   // End of SPD
911
912   
913   //TH1F*  hPrim = (TH1F*)fList->FindObject(fPrimEvents.Data());
914   //  TFile testhit("/home/canute/ALICE/Simulations/TestOfAnalysis2/hitsdist.root");
915   // TH1F* hHitCor = (TH1F*)testhit.Get("selectedHits");
916   // hHitCor->Sumw2();
917   // sumMult->Divide(hHitCor);
918   TH1F*  hPrim = (TH1F*)fList->FindObject(fPrimdNdeta.Data());
919   Float_t nPrimevents = hMCEvents->GetEntries();
920   //std::cout<<hMCEvents->GetEntries()<<std::endl;
921   Float_t xb = hPrim->GetNbinsX();
922   Float_t xr = hPrim->GetXaxis()->GetXmax() - hPrim->GetXaxis()->GetXmin(); 
923   //std::cout<<xb/xr<<std::endl;
924   hPrim->Scale(xb / xr );
925   if(nPrimevents > 0)
926     hPrim->Scale(1/nPrimevents);
927   
928   
929   //  primMult->Rebin(rebin);
930   // primMult->Scale(1/rebin);
931   // hPrim->Rebin(rebin);
932   // hPrim->Scale(1/rebin);
933   if(what != kHits && what != kHitsTrVtx)
934     primMult = hPrim;
935   
936   //  hReldif->Add(hPrim,-1);
937   // hReldif->Divide(hPrim);
938   
939
940
941   /////////////*******************New thing!!!
942   
943   //TH3D* hist3d = (TH3D*)fDataObject->ProjectionXYZ("hist3d");
944   // fDataObject->Sumw2();
945   //TH2F* projeta = (TH2F*)fDataObject->Project3D("yx");
946   
947   //  TProfile2D* projeta = fDataObject->Project3DProfile("yx");
948   // projeta->SetSumw2();
949   //TProfile* etaprofile = projeta->ProfileX("dNdeta_profile");
950   //TProfile* etaprofile = projeta->ProfileX("dNdeta_profile");
951   // TH1* etaprofile = projeta->ProjectionX("dNdeta_profile");
952   //TProfile* etaprofile = fDataObject->ProfileX();
953   /*
954   TCanvas* ctest = new TCanvas("test","test",1200,600);
955   ctest->Divide(3);
956   ctest->cd(1);
957   fDataObject->DrawCopy("box");
958   ctest->cd(2);
959   projeta->DrawCopy("colz");
960   ctest->cd(3);
961   etaprofile->DrawCopy();
962   */
963   //sumMult->Rebin(rebin);
964   TH1F* hReldif = (TH1F*)sumMult->Clone("hReldif");
965   if(rebin != 1) {
966     RebinHistogram(sumMult,rebin);
967     if(!realdata) {
968       RebinHistogram(primMult,rebin);
969       RebinHistogram(hReldif,rebin);
970     }
971   }
972   hReldif->Add(primMult,-1);
973   hReldif->Divide(primMult);
974   TCanvas* c2 = new TCanvas("dN/deta","dN/deta",640,960);
975   c2->Divide(1,2);//,0,0);
976   c2->cd(2);
977   gPad->Divide(1,2);
978   gPad->cd(1);
979   gPad->SetGridy();
980   gPad->SetGridx();
981   hReldif->SetTitle("");
982   hReldif->GetYaxis()->SetTitle("Relative difference");
983   hReldif->GetYaxis()->SetRangeUser(-0.2,0.2);
984   hReldif->SetStats(0000);
985   hReldif->GetXaxis()->SetTitle("#eta");
986   
987   hReldif->DrawCopy();
988   //SPD Rel dif
989   TH1F* hReldifSPD = (TH1F*)hSPDcombi->Clone("hReldifSPD");
990   if(rebin != 1) {
991     RebinHistogram(hSPDcombi,rebin);
992     if(!realdata) {
993       RebinHistogram(hReldifSPD,rebin);
994     }
995   }
996   hReldifSPD->Add(primMult,-1);
997   hReldifSPD->Divide(primMult);
998   hReldifSPD->DrawCopy("same");
999   
1000   TLine* zeroLine = new TLine(-4,0,6,0);
1001   zeroLine->SetLineWidth(2);
1002   zeroLine->SetLineColor(kBlack);
1003   zeroLine->Draw();
1004   hPrim->SetTitle("");
1005   hPrim->SetLabelFont(132,"xyz");
1006   hPrim->SetFillColor(kGray);
1007   primMult->SetFillColor(kBlue);
1008   
1009   
1010   c2->cd(1);
1011   gPad->SetGridy();
1012   gPad->SetGridx();
1013   hPrim->GetYaxis()->SetRangeUser(2,10);
1014   //hPrim->GetYaxis()->SetRangeUser(0,20);
1015   hPrim->SetStats(0000);
1016   hPrim->GetXaxis()->SetTitle("#eta");
1017   sumMult->SetTitle("");
1018   sumMult->SetStats(kFALSE);
1019   sumMult->SetXTitle("#eta");
1020   sumMult->SetYTitle("#frac{dN}{d#eta_{ch}}");
1021   // sumMult->GetYaxis()->SetRangeUser
1022   //sumMult->Scale(1/(Float_t)rebin);
1023   sumMult->DrawCopy("PE");
1024   //Syst errors
1025   TH1F* sumMultSystPos = (TH1F*)sumMult->Clone("systerrorsPos");
1026   TH1F* sumMultSystNeg = (TH1F*)sumMult->Clone("systerrorsNeg");
1027   for(Int_t jj = 1;jj <=sumMultSystPos->GetNbinsX();jj++) {
1028     if(sumMultSystPos->GetBinCenter(jj) < 0) {
1029       sumMultSystPos->SetBinError(jj,0);
1030       sumMultSystPos->SetBinContent(jj,0);
1031       continue;
1032     }
1033       
1034     if(sumMultSystPos->GetBinContent(jj) > 0)
1035       sumMultSystPos->SetBinError(jj,TMath::Sqrt(TMath::Power(sumMultSystPos->GetBinError(jj),2)+TMath::Power(0.1*sumMultSystPos->GetBinContent(jj),2)));
1036   }
1037   for(Int_t jj = 1;jj <=sumMultSystNeg->GetNbinsX();jj++) {
1038     if(sumMultSystNeg->GetBinCenter(jj) > 0) {
1039       sumMultSystNeg->SetBinError(jj,0);
1040       sumMultSystNeg->SetBinContent(jj,0);
1041       continue;
1042     }
1043       
1044     if(sumMultSystNeg->GetBinContent(jj) > 0)
1045       sumMultSystNeg->SetBinError(jj,TMath::Sqrt(TMath::Power(sumMultSystNeg->GetBinError(jj),2)+TMath::Power(0.1*sumMultSystNeg->GetBinContent(jj),2)));
1046   }
1047   
1048   
1049   
1050   sumMultSystPos->SetFillColor(18);
1051   sumMultSystNeg->SetFillColor(18);
1052   //sumMultSystPos->DrawCopy("sameE5");
1053   //sumMultSystNeg->DrawCopy("sameE5");
1054   //End syst errors
1055   
1056   
1057   
1058   if(what != kHits && what != kHitsTrVtx && hPrim->GetEntries())
1059     hPrim->DrawCopy("E3same");
1060   if(what == kHits || what == kHitsTrVtx)
1061     primMult->DrawCopy("E3same");
1062   hPrim->GetXaxis()->SetTitle("#eta");
1063   
1064   TH1F* hMirror = new TH1F("hMirror","mirrored",sumMult->GetNbinsX(),sumMult->GetXaxis()->GetXmin(),sumMult->GetXaxis()->GetXmax());
1065   
1066   for(Int_t i=0.5*sumMult->GetNbinsX(); i<=sumMult->GetNbinsX(); i++) {
1067     Float_t eta   = sumMult->GetBinCenter(i);
1068     Float_t value = sumMult->GetBinContent(i);
1069     Float_t error = sumMult->GetBinError(i);
1070     if(value > 0) {
1071       hMirror->SetBinContent(hMirror->FindBin(-1*eta),value);
1072       hMirror->SetBinError(hMirror->FindBin(-1*eta),error);
1073     }
1074   }
1075   TH1F* hReldifLR = (TH1F*)sumMult->Clone("hReldifLeftRight");
1076   hReldifLR->Add(hMirror,-1);
1077   hReldifLR->Divide(hMirror);
1078   c2->cd(2);
1079   gPad->cd(1);
1080   hReldifLR->GetYaxis()->SetRangeUser(-0.2,0.2);
1081
1082   hReldifLR->SetYTitle("Relative difference left-right");
1083   hReldifLR->SetLabelSize(2*hReldifLR->GetLabelSize("X"),"X"); 
1084   hReldifLR->SetLabelSize(2*hReldifLR->GetLabelSize("Y"),"Y");
1085   hReldifLR->SetTitleSize(2*hReldifLR->GetTitleSize("X"),"X");
1086   hReldifLR->SetTitleSize(2*hReldifLR->GetTitleSize("Y"),"Y");
1087   hReldifLR->SetTitleOffset(0.7*hReldifLR->GetTitleOffset("X"),"X");
1088   hReldifLR->SetTitleOffset(0.5*hReldifLR->GetTitleOffset("Y"),"Y");
1089   
1090   
1091   if(realdata) 
1092     hReldifLR->DrawCopy();
1093   zeroLine->Draw();
1094   c2->cd(1);
1095   hMirror->SetMarkerStyle(25);
1096   
1097   hPrim->SetMarkerStyle(8);
1098   
1099   //graph2->Draw("P");
1100   TH1F* hPythiaMB = 0;
1101   
1102   TFile* fpyt = 0;
1103   
1104   if(!pythiafile.Contains("none"))
1105     fpyt = TFile::Open(pythiafile);
1106   
1107   if(realdata ) {
1108     if(fpyt && !pythiafile.Contains("none")) {
1109       hPythiaMB = (TH1F*)fpyt->Get("hPythiaMB");
1110     }
1111     else
1112       std::cout<<"no pythia for this energy"<<std::endl;
1113   }
1114   if(what != kHits && what != kHitsTrVtx) {
1115     if(hPrim->GetEntries() != 0 && !realdata)
1116       hPrim->DrawCopy("PE3same");
1117     if(realdata) {
1118       //graph->Draw("PEsame");
1119       //graph2->Draw("PEsame");
1120       if(what != kMultNSD) {
1121         graphinel->Draw("PEsame");
1122         graphinel2->Draw("PEsame");
1123         p7742_d1x1y1->Draw("PEsame");
1124       }
1125       else{
1126         graph->Draw("PEsame");
1127         graph2->Draw("PEsame");
1128         p7742_d2x1y1->Draw("PEsame");
1129         p7743_d8x1y1->Draw("PEsame");
1130
1131       }
1132         
1133     }
1134
1135   }
1136   
1137   
1138   sumMult->DrawCopy("PEsame");
1139   
1140   
1141         
1142   hSPDcombi->DrawCopy("same");
1143   if(realdata)
1144     hMirror->DrawCopy("PEsame");
1145     
1146   TLegend* leg = new TLegend(0.3,0.2,0.7,0.45,"");
1147   if(!realdata) {
1148     if(what != kHits && what != kHitsTrVtx)
1149       leg->AddEntry(hPrim,"Primary","pf");
1150     else
1151       leg->AddEntry(primMult,"Hits","pf");
1152   }
1153   
1154   if(what == kMult)
1155     leg->AddEntry(sumMult,"FMD INEL","p");
1156   else if(what == kMultTrVtx)
1157     leg->AddEntry(sumMult,"FMD TrVtx","p");
1158   else if(what == kMultNSD)
1159     leg->AddEntry(sumMult,"FMD NSD","p");
1160   
1161   
1162   if(realdata) {
1163   if(what == kMult)
1164     leg->AddEntry(hMirror,"Mirror FMD INEL","p");
1165   else if(what == kMultTrVtx)
1166     leg->AddEntry(hMirror,"Mirror FMD TrVtx","p");
1167   else if(what == kMultNSD)
1168     leg->AddEntry(hMirror,"Mirror FMD NSD","p");
1169    
1170   if(what == kMultNSD) {  
1171     leg->AddEntry(graph,"UA5 NSD","p");
1172     leg->AddEntry(graph2,"Mirror UA5 NSD","p"); }
1173   else if(what == kMult) {
1174     leg->AddEntry(graphinel,"UA5 INEL","p");
1175     leg->AddEntry(graphinel2,"Mirror UA5 INEL","p");
1176   }
1177     //leg->AddEntry(hPythiaMB,"Pythia MB","l");
1178   leg->AddEntry(hSPDcombi,"HHD SPD clusters","p");
1179   if(what == kMult)
1180     leg->AddEntry(p7742_d1x1y1,"ALICE INEL published","p");
1181   if(what == kMultNSD) {
1182     leg->AddEntry(p7742_d2x1y1,"ALICE NSD published","p");
1183     leg->AddEntry(p7743_d8x1y1, "CMS NSD published","p");
1184   }
1185   }
1186   leg->Draw();
1187   
1188   c2->cd(2);
1189   gPad->cd(2);
1190   TH1F* hRatioMultPythia = 0;
1191   TH1F* hRatioMultUA5 = 0;
1192   TH1F* hRatioMultUA5_rebin5 = new TH1F("hRatioMultUA5_rebin5","hRatioMultUA5_rebin5",sumMult->GetNbinsX(),sumMult->GetXaxis()->GetXmin(),sumMult->GetXaxis()->GetXmax());
1193   
1194   Double_t xval=0, yval=0;
1195   Int_t ipos=0;
1196   for(Int_t bb =hRatioMultUA5_rebin5->FindBin(0); bb<=hRatioMultUA5_rebin5->GetNbinsX();bb++) {
1197     
1198     
1199     
1200     xval=yval=0;
1201     
1202     if(hRatioMultUA5_rebin5->GetBinCenter(bb) >= 0) {
1203       if(what == kMultNSD)
1204         graph->GetPoint(ipos,xval,yval);
1205       else
1206         graphinel->GetPoint(ipos,xval,yval);
1207       
1208       if(yval>0) {
1209         hRatioMultUA5_rebin5->SetBinContent(bb,yval);
1210         hRatioMultUA5_rebin5->SetBinError(bb,graphinel->GetErrorY(ipos));
1211         if(hRatioMultUA5_rebin5->GetBinCenter(bb) < 4) {
1212           hRatioMultUA5_rebin5->SetBinContent(hRatioMultUA5_rebin5->FindBin(-1*hRatioMultUA5_rebin5->GetBinCenter(bb)),yval);
1213           hRatioMultUA5_rebin5->SetBinError(hRatioMultUA5_rebin5->FindBin(-1*hRatioMultUA5_rebin5->GetBinCenter(bb)),(what == kMultNSD ? graph->GetErrorY(ipos) : graph->GetErrorY(ipos)));
1214         
1215         }
1216         ipos++;
1217       }
1218     }
1219   }
1220        
1221       
1222
1223   if(realdata) {
1224     
1225     Float_t ratio = 1;
1226     if(hPythiaMB) {
1227       if(hPythiaMB->GetNbinsX() !=  sumMult->GetNbinsX())
1228         ratio = (Float_t)sumMult->GetNbinsX() / (Float_t)hPythiaMB->GetNbinsX();
1229     }
1230     
1231     hRatioMultPythia = (TH1F*)sumMult->Clone("hRatioMultPythia");
1232     hRatioMultUA5    = (TH1F*)sumMult->Clone("hRatioMultUA5");
1233     if(ratio > 1) {
1234       hRatioMultPythia->Rebin(ratio);
1235       hRatioMultPythia->Scale(1/ratio);
1236     }
1237     if(ratio < 1 && hPythiaMB) {
1238       hPythiaMB->Rebin(1/ratio);
1239       hPythiaMB->Scale(ratio);
1240     }
1241     
1242     if(rebin !=5) {
1243       TGraphErrors* tmp1;
1244       if(what == kMultNSD)
1245       tmp1 = (TGraphErrors*)graph->Clone("UA5tmp");
1246     else
1247       tmp1 = (TGraphErrors*)graphinel->Clone("UA5tmp");
1248     
1249     tmp1->Fit("pol8","Q0","Q0",1.5,6);
1250       TF1* hFit = tmp1->GetFunction("pol8");
1251       for(Int_t ii = hRatioMultUA5->GetNbinsX() / 2; ii<hRatioMultUA5->GetNbinsX();ii++) { 
1252         if(hRatioMultUA5->GetBinContent(ii) > 0) {
1253           Float_t errorratio = hRatioMultUA5->GetBinError(ii) / hRatioMultUA5->GetBinContent(ii);
1254           hRatioMultUA5->SetBinContent(ii,
1255                                        hRatioMultUA5->GetBinContent(ii) / 
1256                                        ((hRatioMultUA5->GetBinCenter(ii) > 4.8) ? hFit->Eval(hRatioMultUA5->GetBinCenter(ii-1)) : hFit->Eval(hRatioMultUA5->GetBinCenter(ii))));
1257           hRatioMultUA5->SetBinError(ii,hRatioMultUA5->GetBinContent(ii) * TMath::Sqrt(0.02*0.02 + TMath::Power(errorratio,2)));
1258           
1259             
1260         }
1261           
1262       }
1263       
1264       TGraphErrors* tmp2;
1265       if(what == kMultNSD)
1266         tmp2 = (TGraphErrors*)graph2->Clone("UA5tmp2");
1267       else 
1268         tmp2 = (TGraphErrors*)graphinel2->Clone("UA5tmp2");
1269
1270       //tmp2->Fit("pol8","Q0","Q0",-3.7,-1.5);
1271       tmp2->Fit("pol7","Q0","Q0",-3.7,0);
1272       hFit = tmp2->GetFunction("pol7");
1273       for(Int_t ii = 0; ii<hRatioMultUA5->GetNbinsX()/2;ii++) { 
1274         if(hRatioMultUA5->GetBinContent(ii) > 0) {
1275           Float_t errorratio = hRatioMultUA5->GetBinError(ii) / hRatioMultUA5->GetBinContent(ii);
1276           hRatioMultUA5->SetBinContent(ii,hRatioMultUA5->GetBinContent(ii) / hFit->Eval(hRatioMultUA5->GetBinCenter(ii))  );
1277           hRatioMultUA5->SetBinError(ii,hRatioMultUA5->GetBinContent(ii) * TMath::Sqrt(0.02*0.02 + TMath::Power(errorratio,2)));
1278           
1279         }
1280         
1281         
1282         graphinel->GetHistogram()->SetStats(kFALSE);
1283         graphinel2->GetHistogram()->SetStats(kFALSE);
1284         
1285       }
1286     }
1287     else hRatioMultUA5->Divide(hRatioMultUA5_rebin5);
1288       //}
1289   
1290       
1291   
1292       gPad->SetGridx();
1293       gPad->SetGridy();
1294     TLine* oneLine = new TLine(-4,1,6,1);
1295     oneLine->SetLineWidth(2);
1296     
1297     hRatioMultUA5->SetMarkerColor(kGreen);
1298     hRatioMultUA5->SetMarkerStyle(20);
1299     hRatioMultUA5->GetYaxis()->SetRangeUser(0.75,1.25);
1300     
1301     hRatioMultUA5->SetYTitle("Ratio FMD/UA5");
1302     hRatioMultUA5->SetLabelSize(2*hRatioMultUA5->GetLabelSize("X"),"X"); 
1303     hRatioMultUA5->SetLabelSize(2*hRatioMultUA5->GetLabelSize("Y"),"Y");
1304     hRatioMultUA5->SetTitleSize(2*hRatioMultUA5->GetTitleSize("X"),"X");
1305     hRatioMultUA5->SetTitleSize(2*hRatioMultUA5->GetTitleSize("Y"),"Y");
1306     hRatioMultUA5->SetTitleOffset(0.7*hRatioMultUA5->GetTitleOffset("X"),"X");
1307     hRatioMultUA5->SetTitleOffset(0.5*hRatioMultUA5->GetTitleOffset("Y"),"Y");
1308     
1309     
1310     
1311     hRatioMultUA5->DrawCopy();
1312     //hRatioMultPythia->DrawCopy("same");
1313     oneLine->Draw("same");
1314   }
1315   
1316   TCanvas* cPaper = new TCanvas("FMD_SPD_UA5","FMD_SPD_UA5",800,600);
1317   cPaper->cd();
1318   
1319   sumMult->DrawCopy();
1320   //sumMultSystPos->DrawCopy("sameE5");
1321   // sumMultSystNeg->DrawCopy("sameE5");
1322   //hTestdNdeta->DrawCopy("same");
1323   //hSPDcombi->DrawCopy("same");
1324   if(what == kMult)
1325     p7742_d1x1y1->Draw("PEsame");
1326   if(what == kMultNSD) {
1327     p7742_d2x1y1->Draw("PEsame");
1328     p7743_d8x1y1->Draw("PEsame");
1329   }
1330
1331   TLegend* leg2 = new TLegend(0.3,0.2,0.7,0.45,"");
1332   if(what == kMult)
1333     leg2->AddEntry(sumMult,"FMD INEL","p");
1334   else if(what == kMultTrVtx)
1335     leg2->AddEntry(sumMult,"FMD TrVtx","p");
1336   else if(what == kMultNSD) 
1337     leg2->AddEntry(sumMult,"FMD NSD","p");
1338   
1339   if(realdata) {
1340     
1341     /* if(what == kMult)
1342     leg2->AddEntry(hMirror,"Mirror FMD INEL","p");
1343   else if(what == kMultTrVtx)
1344     leg2->AddEntry(hMirror,"FMD TrVtx","p");
1345   else if(what == kMultNSD)
1346   leg2->AddEntry(hMirror,"FMD NSD","p");*/
1347    
1348   if(what == kMultNSD) {  
1349     leg2->AddEntry(graph,"UA5 NSD","p");
1350     leg2->AddEntry(graph2,"Mirror UA5 NSD","p"); }
1351   else if(what == kMult) {
1352     leg2->AddEntry(graphinel,"UA5 INEL","p");
1353     leg2->AddEntry(graphinel2,"Mirror UA5 INEL","p");
1354   }
1355     //leg2->AddEntry(hPythiaMB,"Pythia MB","l");
1356     //leg2->AddEntry(hSPDcombi,"HHD SPD clusters","p");
1357   if(what == kMult)
1358     leg2->AddEntry(p7742_d1x1y1,"ALICE INEL published","p");
1359   if(what == kMultNSD) {
1360     leg2->AddEntry(p7742_d2x1y1,"ALICE NSD published","p");
1361     leg2->AddEntry(p7743_d8x1y1,"CMS NSD published","p");
1362   }
1363   }
1364   leg2->Draw();
1365   
1366   
1367   
1368   if(what != kMultNSD) {
1369     graphinel->Draw("PEsame");
1370     graphinel2->Draw("PEsame");
1371   }
1372   else {
1373     graph->Draw("PEsame");
1374     graph2->Draw("PEsame");
1375   }
1376     
1377   sumMult->DrawCopy("PEsame");
1378   c2->cd(1);
1379
1380   c2->Print("fmdana.png");
1381   TString species;
1382   
1383   switch(what) {
1384   case kMult: 
1385     species = "mult";
1386     break;
1387   case kMultTrVtx:
1388     species = "mult_TrVtx";
1389     break;
1390   case kHits:
1391     species = "hits";
1392     break;
1393   case kHitsTrVtx:
1394     species = "hits_TrVtx";
1395     break;
1396   case kMultNSD:
1397     species = "mult_NSD";
1398     break;
1399     
1400   default:
1401     AliWarning("no valid Analysis entry!");
1402     break;
1403   }
1404   TFile fout(Form("fmd_dNdeta_%s.root",species.Data()),"recreate");
1405   if(hRatioMultPythia)
1406     hRatioMultPythia->Write();
1407   hPrim->Write();
1408   graph->Write("UA5nsd");
1409   graph2->Write("UA5mirrornsd");
1410   graphinel->Write("UA5inel");
1411   graphinel2->Write("UA5mirrorinel");
1412   sumMult->Write();
1413   hSPDcombi->Write();
1414   hReldif->Write();
1415   fMultList[what]->Write();
1416   c2->Write();
1417   //fDataObject->Write();
1418   fout.Close();
1419     
1420 }
1421 //_____________________________________________________________________
1422 void AliFMDDndeta::CreateSharingEfficiency(const Char_t* filename, Bool_t store) {
1423   //Generate sharing efficiency
1424   Init(filename);
1425   AliFMDAnaParameters* pars = AliFMDAnaParameters::Instance();
1426   //pars->Init(kTRUE,AliFMDAnaParameters::kBackgroundCorrection);
1427   Int_t nVertexBins = pars->GetNvtxBins();
1428   
1429   SetNames(kHits);
1430   TH1I* hEvents          = (TH1I*)fList->FindObject(fEvents.Data());
1431   TH1I* hPrimEvents      = (TH1I*)fList->FindObject(fPrimEvents.Data());
1432
1433   SetNames(kHitsTrVtx);
1434   TH1I* hEventsTrVtx     = (TH1I*)fList->FindObject(fEvents.Data());
1435   TH1I* hPrimEventsTrVtx = (TH1I*)fList->FindObject(fPrimEvents.Data());
1436   
1437   AliFMDAnaCalibSharingEfficiency* sharEff = new AliFMDAnaCalibSharingEfficiency();
1438   
1439   for(Int_t det = 1; det<=3; det++) {
1440     Int_t maxRing = (det == 1 ? 0 : 1);
1441     for(Int_t ring = 0;ring<=maxRing;ring++) {
1442       Char_t ringChar = (ring == 0 ? 'I' : 'O');
1443       for(Int_t v=0; v< nVertexBins; v++) {
1444         TH1F* hHits = (TH1F*)fList->FindObject(GetAnalysisName(kHits,det, ringChar, v));
1445         TH1F* hHitsTrVtx   = (TH1F*)fList->FindObject(GetAnalysisName(kHitsTrVtx,det, ringChar, v));
1446         TH1F* hMCHits      = (TH1F*)fList->FindObject(GetPrimName(kHits,det, ringChar, v));
1447         TH1F* hMCHitsTrVtx = (TH1F*)fList->FindObject(GetPrimName(kHitsTrVtx,det, ringChar, v));
1448         
1449         TH1F* hCorrection  = (TH1F*)hHits->Clone(Form("hCorrection_FMD%d%c_vtx%d",det, ringChar, v));
1450         TH1F* hCorrectionTrVtx  = (TH1F*)hHitsTrVtx->Clone(Form("hCorrection_FMD%d%c_vtx%d",det, ringChar, v));
1451         if(hEvents->GetBinContent(v+1))
1452           hCorrection->Scale(1./(Float_t)hEvents->GetBinContent(v+1));
1453         if(hPrimEvents->GetBinContent(v+1))
1454           hMCHits->Scale(1./(Float_t)hPrimEvents->GetBinContent(v+1));
1455         if(hEventsTrVtx->GetBinContent(v+1))
1456           hCorrectionTrVtx->Scale(1./(Float_t)hEventsTrVtx->GetBinContent(v+1));
1457         if(hPrimEventsTrVtx->GetBinContent(v+1))
1458           hMCHitsTrVtx->Scale(1./(Float_t)hPrimEventsTrVtx->GetBinContent(v+1));
1459         
1460         hCorrection->Divide(hMCHits);
1461         hCorrectionTrVtx->Divide(hMCHitsTrVtx);
1462         
1463         //sharEff->SetSharingEff(det,ringChar,v,hCorrection);
1464         sharEff->SetSharingEff(det,ringChar,v,hCorrection);
1465         sharEff->SetSharingEffTrVtx(det,ringChar,v,hCorrectionTrVtx);
1466         
1467         //      std::cout<<hHits<<"  "<<hHitsTrVtx<<"   "<<hPrim<<"    "<<hPrimTrVtx<<std::endl;
1468
1469         }
1470       }
1471     }
1472
1473   if(store) {
1474     TFile fsharing(pars->GetPath(pars->GetSharingEffID()),"RECREATE");
1475     sharEff->Write(AliFMDAnaParameters::GetSharingEffID());
1476     fsharing.Close();
1477   }
1478   
1479 }
1480 //_____________________________________________________________________
1481 void AliFMDDndeta::RebinHistogram(TH1F* hist, Int_t rebin) {
1482   //Clever rebinner
1483   
1484   if(hist->GetNbinsX()%rebin) {
1485     std::cout<<"cannot rebin "<<hist->GetNbinsX()%rebin<< "is not zero"<<std::endl;
1486     return;
1487
1488   }
1489   TH1F* cloneHist = (TH1F*)hist->Clone();
1490   cloneHist->Rebin(rebin);
1491   
1492   Int_t nBinsNew = hist->GetNbinsX() / rebin;
1493   
1494   for(Int_t i =1;i<=nBinsNew; i++) {
1495     Float_t bincontent = 0;
1496     Float_t sumofw = 0;
1497     Float_t sumformean = 0;
1498     Int_t   nBins = 0;
1499     for(Int_t j=1; j<=rebin;j++) {
1500       if(hist->GetBinContent((i-1)*rebin + j) > 0) {
1501         bincontent = bincontent + hist->GetBinContent((i-1)*rebin + j);
1502         nBins++;
1503         Float_t weight = 1 / TMath::Power(hist->GetBinError((i-1)*rebin + j),2);
1504         sumofw = sumofw + weight;
1505         sumformean = sumformean + weight*hist->GetBinContent((i-1)*rebin + j);
1506         
1507       }
1508     }
1509     
1510     if(bincontent > 0 ) {
1511       cloneHist->SetBinContent(i,sumformean / sumofw);
1512       cloneHist->SetBinError(i,TMath::Sqrt(sumformean));//TMath::Sqrt(1/sumofw));
1513     }
1514   }
1515   hist->Rebin(rebin);
1516   for(Int_t i =1;i<=nBinsNew; i++) {
1517     hist->SetBinContent(i,cloneHist->GetBinContent(i));
1518   }
1519   
1520   
1521 }
1522 //_____________________________________________________________________
1523 //
1524 // EOF
1525 //