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