Coverity fixes
[u/mrichter/AliRoot.git] / PWG2 / EVCHAR / AliCentralityGlauberFit.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /*   Origin: Alberica Toia, CERN, Alberica.Toia@cern.ch                   */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 //                                                                           //
20 //  class to determine centrality percentiles from 1D distributions          // 
21 //                                                                           //
22 ///////////////////////////////////////////////////////////////////////////////
23
24 #include <TNamed.h>
25 #include <TH1D.h>
26 #include <TString.h>
27 #include <TFile.h>
28 #include <TMath.h>
29 #include <TROOT.h>
30 #include <TH2F.h>
31 #include <TF1.h>
32 #include <TNtuple.h>
33 #include <TStyle.h>
34 #include <TGraphErrors.h>
35 #include <vector>
36 #include <TMinuit.h>
37 #include <TStopwatch.h>
38 #include <TRandom.h>
39 #include "AliCentralityGlauberFit.h"
40 #include "AliLog.h"
41
42 ClassImp(AliCentralityGlauberFit)  
43  
44 //______________________________________________________________________________
45 AliCentralityGlauberFit::AliCentralityGlauberFit(const char *filename) : 
46   fNmu(0),    
47   fMulow(0),  
48   fMuhigh(0), 
49   fNk(0),     
50   fKlow(0),   
51   fKhigh(0),  
52   fNalpha(0), 
53   fAlphalow(0),   
54   fAlphahigh(0),  
55   fRebinFactor(0),
56   fScalemin(0),    
57   fMultmin(0),    
58   fMultmax(0),    
59   fGlauntuple(0), 
60   fNpart(0),       
61   fNcoll(0),       
62   fB(0),           
63   fTaa(0),         
64   fEffi(0),        
65   fhEffi(0),      
66   fTempHist(0),   
67   fGlauberHist(0), 
68   fFastFit(0),      
69   fAncestor(2),      
70   fNBD(0),         
71   fUseChi2(kTRUE),      
72   fUseAverage(kFALSE),
73   fhAncestor(0),
74   fNevents(100000),
75   fNtrials(100),
76   fInrootfilename(0),
77   fInntuplename(0),
78   fOutrootfilename(0),
79   fOutntuplename(0),
80   fAncfilename("ancestor_hists.root"),
81   fHistnames()
82 {
83   // Standard constructor.
84   TFile *f = 0;
85   if (filename) {  // glauber file
86     f = TFile::Open(filename);
87     fGlauntuple = (TNtuple*) f->Get("nt_Pb_Pb");
88     fGlauntuple->SetBranchAddress("Npart",&fNpart);
89     fGlauntuple->SetBranchAddress("Ncoll",&fNcoll);
90     fGlauntuple->SetBranchAddress("B",&fB);
91     fGlauntuple->SetBranchAddress("tAA",&fTaa);
92   }
93
94   fNBD = new TF1 ("fNBD", AliCentralityGlauberFit::NBDFunc, 0, 100,2);
95 }
96
97 //--------------------------------------------------------------------------------------------------
98 void AliCentralityGlauberFit::SetRangeToFit(Double_t fmultmin, Double_t fmultmax)
99 {
100   // Set fit range.
101
102   fMultmin=fmultmin;
103   fMultmax=fmultmax;
104 }
105
106 //--------------------------------------------------------------------------------------------------
107 void AliCentralityGlauberFit::SetRangeToScale(Double_t fscalemin)
108 {
109   // Set range where to scale simulated histo to data
110
111   fScalemin=fscalemin;
112 }
113
114 //--------------------------------------------------------------------------------------------------
115 void AliCentralityGlauberFit::SetGlauberParam(
116                                               Int_t   Nmu,
117                                               Double_t mulow,
118                                               Double_t muhigh,
119                                               Int_t   Nk,
120                                               Double_t klow,
121                                               Double_t khigh,
122                                               Int_t   Nalpha,
123                                               Double_t alphalow,
124                                               Double_t alphahigh)
125 {
126   // Set Glauber parameters.
127
128   fNmu=Nmu;
129   fMulow=mulow;
130   fMuhigh=muhigh;
131   fNk=Nk;
132   fKlow=klow;
133   fKhigh=khigh;
134   fNalpha=Nalpha;
135   fAlphalow=alphalow;
136   fAlphahigh=alphahigh;
137 }
138
139 //--------------------------------------------------------------------------------------------------
140 void AliCentralityGlauberFit::MakeFits() 
141 {
142   // Make fits.
143
144   TH1F *hDATA;
145   TH1F *thistGlau; 
146   TFile *inrootfile;
147   TFile *outrootfile;  
148   //FILE* fTxt = fopen ("parameters.txt","w");
149   
150   // open inrootfile, outrootfile
151   std::cout << "input file "  << fInrootfilename  << std::endl;
152   std::cout << "output file " << fOutrootfilename << std::endl;
153   inrootfile  = new TFile(fInrootfilename,"OPEN");
154   outrootfile = new TFile(fOutrootfilename,"RECREATE");
155   
156   // loop over all distribution names  
157   std::vector<TString>::const_iterator hni;
158   for(hni=fHistnames.begin(); hni!=fHistnames.end(); hni++) {
159     hDATA  = (TH1F*) (inrootfile->Get(*hni)); 
160     if (!hDATA) {
161       TList *list  = (TList*) (inrootfile->Get("CentralityStat")); 
162       hDATA  = (TH1F*) (list->FindObject(*hni));
163     } 
164     hDATA->Rebin(fRebinFactor);
165     //TH1F *hGLAU = new TH1F("hGLAU","hGLAU",hDATA->GetNbinsX(),0,hDATA->GetNbinsX()*hDATA->GetBinWidth(1));
166     Double_t chi2min = 9999999.0;
167     Double_t alpha_min=-1;
168     Double_t mu_min=-1;
169     Double_t k_min=-1;
170     Double_t alpha, mu, k, chi2;   
171
172     for (Int_t nalpha=0;nalpha<fNalpha;nalpha++) {
173       alpha = fAlphalow   + ((Double_t) nalpha ) * (fAlphahigh  - fAlphalow ) / fNalpha;
174
175       mu=0.0;
176       for (Int_t nmu=0; nmu<fNmu; nmu++) {
177         mu = fMulow   + ((Double_t) nmu ) * (fMuhigh  - fMulow ) / fNmu;
178
179         for (Int_t nk=0; nk<fNk; nk++) {
180           k = fKlow  + ((Double_t) nk ) * (fKhigh  - fKlow ) / fNk;
181                     
182           thistGlau = GlauberHisto(mu,k,alpha,hDATA,kFALSE);
183           chi2 = CalculateChi2(hDATA,thistGlau);
184           //fprintf(fTxt, "%3.3f %3.3f %3.3f %3.3f\n", (float) mu, (float) k, (float) alpha, chi2);
185           if ( chi2 < chi2min ) {
186             chi2min = chi2;
187             alpha_min=alpha;
188             mu_min=mu;
189             k_min=k;
190           }
191         }
192       }
193     }
194     
195     thistGlau = GlauberHisto(mu_min,k_min,alpha_min,hDATA,kTRUE);
196     
197     TH1F * hGLAU = 0x0;
198     hGLAU = (TH1F *) thistGlau->Clone("hGLAU");
199     hGLAU->SetName( ((TString)hDATA->GetName()).Append(Form("_GLAU")));
200     hGLAU->SetTitle( ((TString)hDATA->GetName()).Append(Form("_GLAU_%.3f_%.3f_%.3f",
201                                                              mu_min,k_min,alpha_min)));
202
203     Double_t mcintegral = hGLAU->Integral(hGLAU->FindBin(fScalemin),hGLAU->GetNbinsX());
204     Double_t scale = (hDATA->Integral(hDATA->FindBin(fScalemin),hDATA->GetNbinsX())/mcintegral);
205     hGLAU->Scale(scale);
206
207     fhEffi = GetTriggerEfficiencyFunction(hDATA, hGLAU);
208     SaveHisto(hDATA,hGLAU,fhEffi,outrootfile);
209     //fclose (fTxt);
210
211     std::cout << "chi2 min is " << chi2min << std::endl;
212     std::cout << "fitted " << hGLAU->Integral(hGLAU->FindBin(fMultmin),
213                                               hGLAU->FindBin(fMultmax))/hGLAU->Integral() 
214               << " of the total cross section" << std::endl; 
215     fTempHist=hDATA;
216     fGlauberHist=hGLAU;
217   }
218
219   // close inrootfile, outrootfile
220   inrootfile->Close();
221   outrootfile->Close();
222 }
223
224 //--------------------------------------------------------------------------------------------------
225 void AliCentralityGlauberFit::MakeFitsMinuitNBD(Double_t alpha, Double_t mu, Double_t k) 
226 {
227   // Make fits using Minuit.
228
229   if (alpha<0) 
230     alpha = fAlphalow;
231   if (mu<0)
232     mu = fMulow;
233   if (k<0)
234     k = fKlow;
235   printf("Calling Minuit with starting values: %f %f %f\n", alpha, mu, k);
236
237   TH1F *hDATA;
238   TH1F *thistGlau; 
239   TFile *inrootfile;
240   TFile *outrootfile;
241   
242   // open inrootfile, outrootfile
243   std::cout << "input file "  << fInrootfilename  << std::endl;
244   std::cout << "output file " << fOutrootfilename << std::endl;
245   inrootfile  = new TFile(fInrootfilename,"OPEN");
246   outrootfile = new TFile(fOutrootfilename,"RECREATE");
247
248   // loop over all distribution names  
249   std::vector<TString>::const_iterator hni;
250   for(hni=fHistnames.begin(); hni!=fHistnames.end(); hni++) {
251     hDATA  = (TH1F*) (inrootfile->Get(*hni)); 
252     hDATA->Rebin(fRebinFactor);
253     fTempHist=hDATA;
254     TH1F *hGLAU = new TH1F("hGLAU","hGLAU",hDATA->GetNbinsX(),0,hDATA->GetNbinsX()*hDATA->GetBinWidth(1));
255     hGLAU->Sumw2();
256
257     // Minimize here
258     if(gMinuit) delete gMinuit;
259     new TMinuit(3);
260     gMinuit->mncler();
261     gMinuit->SetFCN(AliCentralityGlauberFit::MinuitFcnNBD);
262     gMinuit->SetObjectFit(this); 
263
264     Double_t arglist[2]={0};
265     Int_t ierflg;
266     if (fUseChi2) arglist[0] = 1; // should be 1 if you want to minimize chi2
267     else arglist[0] = 0.5;        // should be 0.5 for ll
268     gMinuit->mnexcm("SET ERR",arglist, 1, ierflg);
269     gMinuit->mnparm(0,"alpha", alpha,  (fAlphahigh-fAlphalow)/fNalpha,  fAlphalow, fAlphahigh, ierflg);
270     gMinuit->mnparm(1,"mu"   , mu,     (fMuhigh-fMulow)/fNmu, fMulow, fMuhigh, ierflg);
271     gMinuit->mnparm(2,"k"    , k,      (fKhigh-fKlow)/fNk, fKlow, fKhigh, ierflg);
272
273     // Call migrad
274     arglist[0] = 100; // max calls
275     arglist[1] = 0.1; // tolerance    
276     gMinuit->mnexcm("SIMPLEX",arglist, 2, ierflg);
277     arglist[0] = 1000; // max calls
278     arglist[1] = 0.1;  // tolerance    
279     gMinuit->mnexcm("MIGrad",arglist, 2, ierflg);
280     //gMinuit->mnexcm("IMPROVE",arglist, 0, ierflg);
281
282     if (ierflg != 0) {
283       AliWarning("Abnormal termination of minimization.");
284     }
285     
286     // ______________________ Get chi2 and Fit Status __________
287     Double_t amin,edm,errdef;
288     Int_t    nvpar, nparx, icstat;
289     
290     gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
291     Double_t chi2min = amin;
292     std::cout << "Fit status " << icstat << std::endl;
293
294     Double_t alpha_min, mu_min, k_min;
295     Double_t alpha_mine, mu_mine, k_mine;
296     gMinuit->GetParameter(0, alpha_min , alpha_mine );
297     gMinuit->GetParameter(1, mu_min    , mu_mine    );
298     gMinuit->GetParameter(2, k_min     , k_mine     );
299
300     // print some infos
301     std::cout << "chi2 min is " << chi2min << ", " << alpha_min << ", "<< mu_min<< ", "  
302               << k_min << std::endl;
303
304     thistGlau = GlauberHisto(mu_min,k_min,alpha_min,hDATA,kTRUE);
305
306     hGLAU = (TH1F *) thistGlau->Clone("hGLAU");
307     hGLAU->SetName( ((TString)hDATA->GetName()).Append(Form("_GLAU")));
308     hGLAU->SetTitle( ((TString)hDATA->GetName()).Append(Form("_GLAU_%.3f_%.3f_%.3f",
309                                                              mu_min,k_min,alpha_min)));
310
311     
312     std::cout << "fitted " << hGLAU->Integral(hGLAU->FindBin(fMultmin),
313                                               hGLAU->FindBin(fMultmax))/hGLAU->Integral() 
314               << " of the total cross section" << std::endl; 
315
316     Double_t mcintegral = hGLAU->Integral(hGLAU->FindBin(fScalemin),hGLAU->GetNbinsX());
317     Double_t scale = (hDATA->Integral(hDATA->FindBin(fScalemin),hDATA->GetNbinsX())/mcintegral);
318     hGLAU->Scale(scale);
319
320     std::cout << "Chi2 final " << CalculateChi2(hDATA,hGLAU) << std::endl;
321
322     fhEffi = GetTriggerEfficiencyFunction(hDATA, hGLAU);
323     SaveHisto(hDATA,hGLAU,fhEffi,outrootfile);
324     fGlauberHist=hGLAU;
325   }
326   // close inrootfile, outrootfile
327   inrootfile->Close();
328   outrootfile->Close();
329 }
330
331 //--------------------------------------------------------------------------------------------------
332 TH1F *AliCentralityGlauberFit::GlauberHisto(Double_t mu, Double_t k, Double_t alpha, 
333                                             TH1F *hDATA, Bool_t save) 
334 {
335   // Get Glauber histogram.
336   static TH1F *h1 = (TH1F*)hDATA->Clone();
337   h1->Reset();
338   h1->SetName(Form("fit_%.3f_%.3f_%.3f",mu,k,alpha));
339
340   if (fUseAverage) {
341     fhAncestor = MakeAncestor(alpha);
342     for (Int_t np=1; np<=fhAncestor->GetNbinsX(); ++np) {  
343       Double_t nanc = fhAncestor->GetBinCenter(np);
344       Double_t weights = fhAncestor->GetBinContent(np);
345
346       if (weights <= 0) continue;
347       Int_t trials = (Int_t) (20 * nanc * (int) mu);
348       if (trials <=0) continue;
349       for (Int_t j=0; j<trials; j++) {
350         double nbdvalue = NBD(j, mu * nanc, k * nanc);
351         h1->Fill((double) j, nbdvalue * weights);
352       }
353     }
354     return h1;
355   }
356
357   TH1F *hSample = fFastFit != 2 ? NBDhist(mu,k) : 0;
358
359   TFile *outFile = NULL;
360   TNtuple *ntuple = NULL;
361  
362   if (save) {
363     outFile = new TFile(fOutntuplename,"RECREATE");
364     ntuple = new TNtuple("gnt", "Glauber ntuple", "Npart:Ncoll:B:tAA:ntot");
365   } 
366
367   Int_t nents = 0;
368   if (fGlauntuple)
369     nents = fGlauntuple->GetEntries(); 
370
371   for (Int_t i=0;i<fNevents;++i) {
372     if (fGlauntuple)
373       fGlauntuple->GetEntry(i % nents);
374     else {
375       fNpart = 2;
376       fNcoll = 1;
377     }
378     Int_t n=0;
379     //if (fAncestor == 1)      n = (Int_t)(TMath::Power(fNpart,alpha));
380     if (fAncestor == 1)      n = (Int_t)(TMath::Power(fNcoll,alpha));
381     else if (fAncestor == 2) n = (Int_t)(alpha * fNpart + (1-alpha) * fNcoll);
382     else if (fAncestor == 3) n = (Int_t)((1-alpha) * fNpart/2 + alpha * fNcoll);
383
384     Int_t ntot=0;
385     if (fFastFit == 1) {
386       ntot = (Int_t)(n*hSample->GetRandom()); // NBD
387     }
388     else if (fFastFit == 2) {
389       Double_t sigma = k*TMath::Sqrt(n*mu);    
390       ntot = (Int_t)(gRandom->Gaus(n*mu,sigma)); // Gaussian
391     }
392     else {
393       for(Int_t j = 0; j<(Int_t)n; ++j) 
394         ntot += (Int_t)hSample->GetRandom();
395     }
396     h1->Fill(ntot);    
397
398     if (save) 
399       ntuple->Fill(fNpart,fNcoll,fB,fTaa,ntot);
400   }
401
402   if (save) {
403     ntuple->Write();
404     outFile->Close();
405   }
406   
407   if (fFastFit != 2) delete hSample;
408   return h1;
409 }
410
411 //--------------------------------------------------------------------------------------------------
412 Double_t AliCentralityGlauberFit::CalculateChi2(TH1F *hDATA, TH1F *thistGlau) 
413 {
414   // Note that for different values of neff the mc distribution is identical, just re-scaled below
415   // normalize the two histogram, scale MC up but leave data alone - that preserves error bars !!!!
416   
417   Int_t lowchibin =   hDATA->FindBin(fMultmin);
418   Int_t highchibin =  hDATA->FindBin(fMultmax);
419
420   Double_t mcintegral = thistGlau->Integral(1,thistGlau->GetNbinsX());
421   Double_t scale = (hDATA->Integral(1,hDATA->GetNbinsX())/mcintegral);
422   thistGlau->Scale(scale);
423
424   // calculate the chi2 between MC and real data over some range ????
425   if (fUseChi2) {
426     Double_t chi2 = 0.0;
427     for (Int_t i=lowchibin; i<=highchibin; i++) {
428       if (hDATA->GetBinContent(i) < 1.0) continue;
429       Double_t diff = TMath::Power((thistGlau->GetBinContent(i) - hDATA->GetBinContent(i)),2);
430       diff = diff / (TMath::Power(hDATA->GetBinError(i),2) + TMath::Power(thistGlau->GetBinError(i),2)); 
431       chi2 += diff;
432     }
433     chi2 = chi2 / (highchibin - lowchibin + 1);
434     return chi2;
435   }
436   // "-2 log likelihood ratio(mu;n) = 2[(mu - n) + n ln(n/mu)]"
437   else {
438     std::cout << "LL" << std::endl;
439
440     Double_t ll = 0.0;
441     for (Int_t i=lowchibin; i<=highchibin; i++) {
442       Double_t data = hDATA    ->GetBinContent(i);
443       Double_t mc   = thistGlau->GetBinContent(i);
444       Int_t idata = TMath::Nint(data);
445       if (mc < 1.e-9) mc = 1.e-9;
446       Double_t fsub = - mc + idata * TMath::Log(mc);
447       Double_t fobs = 0;
448       if (idata > 0) {
449         for(Int_t istep = 0; istep < idata; istep++){
450           if (istep > 1) 
451             fobs += TMath::Log(istep);
452         }
453       }
454       fsub -= fobs;
455       ll -=  fsub ;
456     }
457     return 2*ll;
458   }
459 }
460
461 //--------------------------------------------------------------------------------------------------
462 TH1F *AliCentralityGlauberFit::GetTriggerEfficiencyFunction(TH1F *hist1, TH1F *hist2) 
463 {
464   // Get efficiency.
465
466   fhEffi = (TH1F*)hist1->Clone("heffi");
467   fhEffi->Divide(hist2);
468   return fhEffi;
469 }
470
471 //--------------------------------------------------------------------------------------------------
472 Double_t AliCentralityGlauberFit::GetTriggerEfficiencyIntegral(TH1F *hist1, TH1F *hist2) 
473 {
474   // Get eff integral.
475
476   fEffi = hist1->Integral()/hist2->Integral();
477   return fEffi;
478 }
479
480 //--------------------------------------------------------------------------------------------------
481 void AliCentralityGlauberFit::SaveHisto(TH1F *hist1, TH1F *hist2, TH1F *heffi, TFile *outrootfile) 
482 {
483   // Save histograms.
484
485   outrootfile->cd();
486   hist1->Write();
487   hist2->Write();
488   heffi->Write();
489 }
490
491 //--------------------------------------------------------------------------------------------------
492 Double_t AliCentralityGlauberFit::NBD(Int_t n, Double_t mu, Double_t k) const
493 {
494   // Compute NBD.
495   double F;
496   double f;
497
498   if (n+k > 100.0) {
499     // log method for handling large numbers
500     F  = TMath::LnGamma(n + k)- TMath::LnGamma(n + 1.)- TMath::LnGamma(k);
501     f  = n * TMath::Log(mu/k) - (n + k) * TMath::Log(1.0 + mu/k);
502     F = F+f;
503     F = TMath::Exp(F);
504   } else {
505     F  = TMath::Gamma(n + k) / ( TMath::Gamma(n + 1.) * TMath::Gamma(k) );
506     f  = n * TMath::Log(mu/k) - (n + k) * TMath::Log(1.0 + mu/k);
507     f  = TMath::Exp(f);
508     F *= f;
509   }
510
511   return F;
512
513 }
514
515 //--------------------------------------------------------------------------------------------------
516 TH1F *AliCentralityGlauberFit::NBDhist(Double_t mu, Double_t k)
517 {
518   // Interface for TH1F.
519
520   TH1F *h = new TH1F("htemp","",100,-0.5,299.5);
521   h->SetName(Form("nbd_%f_%f",mu,k));
522   h->SetDirectory(0);
523   for (Int_t i=0;i<300;++i) {
524     Double_t val = NBD(i,mu,k);
525     if (val>1e-20) h->Fill(i,val);
526   }
527   return h;
528 }
529
530 //--------------------------------------------------------------------------------------------------
531 void AliCentralityGlauberFit::MinuitFcnNBD(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
532 {
533   // FCN for minuit.
534
535   Double_t alpha = par[0];
536   Double_t mu    = par[1];
537   Double_t k     = par[2];
538
539   if (0) { //avoid warning
540     gin=gin;
541     npar=npar;
542     iflag=iflag;
543   }
544   AliCentralityGlauberFit * obj = (AliCentralityGlauberFit *) gMinuit->GetObjectFit();
545   TH1F * thistGlau = obj->GlauberHisto(mu,k,alpha,obj->GetTempHist(),kFALSE);
546   f = obj->CalculateChi2(obj->GetTempHist(),thistGlau);
547   Printf("Minuit step: chi2=%f, alpha=%f,  mu=%f, k=%f\n",f,alpha,mu,k);
548 }
549
550 //--------------------------------------------------------------------------------------------------
551 Double_t AliCentralityGlauberFit::NBDFunc(const Double_t *x, const Double_t *par) 
552 {
553   // TF1  interface.
554
555   Double_t mu = par[0];
556   Double_t k  = par[1];
557   Double_t n  = x[0];
558   Double_t ret = exp( TMath::LnGamma(n+k) - TMath::LnGamma(k) - TMath::LnGamma(n+1) ) * 
559     TMath::Power(mu/(mu+k),n) * TMath::Power(1-mu/(mu+k),k);
560   return ret;
561 }
562
563 //--------------------------------------------------------------------------------------------------
564 TH1F *AliCentralityGlauberFit::MakeAncestor(Double_t alpha)
565 {
566   // Make ancestor histogram.
567
568   TString hname(Form("fhAncestor_%.3f",alpha));
569   if (fhAncestor) {
570     if (hname.CompareTo(fhAncestor->GetName())==0)
571       return fhAncestor;
572   }
573
574   delete fhAncestor;
575   fhAncestor = 0;
576   TFile *ancfile = TFile::Open(fAncfilename,"read");
577   if (ancfile && ancfile->IsOpen()) {
578     fhAncestor = dynamic_cast<TH1F*>(ancfile->Get(hname));
579     if (fhAncestor) {
580       fhAncestor->SetDirectory(0);
581       delete ancfile;
582       return fhAncestor;
583     }
584   }
585   delete ancfile;
586   
587   fhAncestor = new TH1F(hname,hname,3000,0,3000);
588   fhAncestor->SetDirectory(0);
589   Int_t nents = fGlauntuple->GetEntries(); 
590   for (Int_t i=0;i<nents;++i) {
591     fGlauntuple->GetEntry(i % nents);
592     Int_t n=0;
593     if (fAncestor == 1)    n = (Int_t) (TMath::Power(fNpart,alpha));
594     //if (fAncestor == 1)      n = (Int_t) (TMath::Power(fNcoll,alpha));
595     else if (fAncestor == 2) n = (Int_t) (alpha * fNpart + (1-alpha) * fNcoll);
596     else if (fAncestor == 3) n = (Int_t) ((1-alpha) * fNpart/2 + alpha * fNcoll);
597     fhAncestor->Fill(n);
598   }
599
600   ancfile = TFile::Open(fAncfilename,"update");
601   if (ancfile && ancfile->IsOpen()) {
602     fhAncestor->Write();
603   }
604   delete ancfile;
605   return fhAncestor;
606 }