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