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