PWG4/macros/runJetSpectrumUnfolding.C

   1 //
   2 // script with functions to use AliJetSpectrumUnfolding class
   3 //
   4
   5
   6 void loadlibs(){
   7   // load all the needed libs to run wit root only
   8
   9   gSystem->Load("libTree");
  10   gSystem->Load("libPhysics");
  11   gSystem->Load("libHist");
  12   gSystem->Load("libVMC");
  13   gSystem->Load("libSTEERBase");
  14   gSystem->Load("libAOD");
  15   gSystem->Load("libESD");
  16   gSystem->Load("libANALYSIS");
  17   gSystem->Load("libANALYSISalice");
  18   gSystem->Load("libJETAN");
  19   gSystem->Load("libPWG4JetTasks");
  20
  21
  22 }
  23
  24
  25
  26 void draw(const char* fileName = "unfolded_pwg4spec.root", const char* folder = "unfolding", Bool_t proj = kTRUE)
  27 {
  28
  29   loadlibs();
  30
  31   AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder);
  32
  33   TFile::Open(fileName);
  34   jetSpec->LoadHistograms();
  35
  36
  37   if (proj)
  38   {
  39     canvas1 = new TCanvas("Response Map Projection", "Response Map Projection", 500, 500);
  40     canvas1->Divide(2);
  41
  42     Int_t style = 1;
  43     const Int_t NRGBs = 5;
  44     const Int_t NCont = 500;
  45
  46     Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
  47     Double_t red[NRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51 };
  48     Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
  49     Double_t blue[NRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00 };
  50     TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
  51     gStyle->SetNumberContours(NCont);
  52     gStyle->SetPalette(style);
  53
  54     canvas1->cd(1);
  55     gPad->SetLogz();
  56     h2 = (jetSpec->GetCorrelation())->Projection(2,3);
  57     h2->SetXTitle("z^{lp}_{gen}");
  58     h2->SetYTitle("z^{lp}_{rec}");
  59     h2->DrawCopy("colz");
  60
  61     canvas1->cd(2);
  62     gPad->SetLogz();
  63     h3 = (jetSpec->GetCorrelation())->Projection(0,1);
  64     h3->SetXTitle("E^{jet}_{gen} [GeV]");
  65     h3->SetYTitle("E^{jet}_{rec} [GeV]");
  66     h3->DrawCopy("colz");
  67   }
  68   jetSpec->DrawComparison("Draw_unfolded_pwg4spec", jetSpec->GetGenSpectrum());
  69
  70   return;
  71 }
  72
  73 //________________________________________________________________________________________________________________
  74 void unfold(const char* fileNameGen = "gen_pwg4spec.root", const char* folder = "unfolding", const char* fileNameRec = "rec_pwg4spec.root", const char* fileNameUnf = "unfolded_pwg4spec.root")
  75 {
  76   // function to load jet spectra from the output file of the task AliAnalysisTaskJetSpectrum
  77   // to do the unfolding
  78
  79   loadlibs();
  80
  81   AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder);
  82
  83   TFile::Open(fileNameGen);
  84   jetSpec->LoadHistograms();
  85
  86   TFile::Open(fileNameRec);
  87   TH2F* hist = (TH2F*) gFile->Get("unfolding/fRecSpectrum");
  88   jetSpec->SetRecSpectrum(hist);
  89   hist = (TH2F*) gFile->Get("unfolding/fGenSpectrum");
  90   if(hist->GetEntries()>0)jetSpec->SetGenSpectrum(hist);
  91
  92   jetSpec->ApplyBayesianMethod(0.3, 20, 0, 0);
  93   // last parameter = calculateErrors  <- this method to calculate the errors takes a lot of time
  94
  95   TFile* file = TFile::Open(fileNameUnf,"RECREATE");
  96   jetSpec->SaveHistograms();
  97   file->Close();
  98 }
  99
 100 //___________________________________________________________________________
 101 void buildSpectra(Int_t caseNo, const char* inFile, const char* outFile)
 102 {
 103   // function to test 2D Bayesian unfolding with other spectra
 104   // build from a function
 105
 106   loadlibs();
 107
 108
 109
 110   AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding("unfolding", "unfolding");
 111
 112   TFile::Open(inFile);
 113   jetSpec->LoadHistograms("unfolding");
 114
 115   TCanvas *c1 = new TCanvas();TH2 *h2 = (jetSpec->GetCorrelation())->Projection(0,3);
 116   h2->DrawCopy("colz");
 117   c1->Update();
 118   if(getchar()=='q')return;
 119
 120
 121   switch (caseNo)
 122   {
 123     case 1: func = new TF2("func", "501-x-y"); break;
 124     case 2: func = new TF2("func", "1000 * 1/(y+x+1)"); break;
 125     case 3: func = new TF2("func", "1./(y*pow(x,5.7))"); break;
 126     case 4: func = new TF2("func", "exp(-0.1*x - 0.1*y)"); break;
 127     case 5: func = new TF2("func", "x*x + (y**3)/100."); break;
 128     case 6: func = new TF2("func", "1000*y*exp(-0.1*x)"); break;
 129     case 7: func = new TF2("func", "exp(-((x-100.)/(0.3*x))**2 - ((y-0.6)/(0.8*y))**2)"); break;
 130     default: return;
 131   }
 132
 133   //new TCanvas; func->Draw();
 134
 135   jetSpec->SetGenRecFromFunc(func);
 136
 137   TFile* file = TFile::Open(outFile,"RECREATE");
 138   jetSpec->SaveHistograms();
 139
 140   h2 = (jetSpec->GetCorrelation())->Projection(0,3);
 141   h2->DrawCopy("colz");
 142   c1->Update();
 143   file->Close();
 144
 145   //new TCanvas; jetSpec->GetRecSpectrum()->DrawCopy("colz");
 146   //new TCanvas; jetSpec->GetGenSpectrum()->DrawCopy("colz");
 147 }
 148
 149 //_____________________________________________________________________________________________
 150 void buildResponseMap(const char* fileName = "responseMap.root")
 151 {
 152   // function to build a Response Map with a gaussian distribution
 153   loadlibs();
 154
 155   AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding("unfolding", "unfolding");
 156
 157   TF2* func = new TF2("func", "exp(-((x-[0])/[1])**2 - ((y-[2])/[3])**2)");
 158
 159   bool bPerfect = false;
 160
 161   Double_t var[4];
 162   Float_t sigmax, sigmay;
 163   for (Int_t tx=1; tx<=jetSpec->GetCorrelation()->GetAxis(0)->GetNbins(); tx++)
 164     for (Int_t ty=1; ty<=jetSpec->GetCorrelation()->GetAxis(2)->GetNbins(); ty++)
 165     {
 166       var[0] = jetSpec->GetCorrelation()->GetAxis(0)->GetBinCenter(tx);
 167       var[2] = jetSpec->GetCorrelation()->GetAxis(2)->GetBinCenter(ty);
 168       sigmax = 0.2*var[0];
 169       sigmay = 0.2*var[2];
 170       func->SetParameters(var[0],sigmax,var[2],sigmay);
 171       for (Int_t mx=1; mx<=jetSpec->GetCorrelation()->GetAxis(1)->GetNbins(); mx++)
 172         for (Int_t my=1; my<=jetSpec->GetCorrelation()->GetAxis(3)->GetNbins(); my++)
 173         {
 174           var[1] = jetSpec->GetCorrelation()->GetAxis(1)->GetBinCenter(mx);
 175           var[3] = jetSpec->GetCorrelation()->GetAxis(3)->GetBinCenter(my);
 176
 177
 178           if(bPerfect){
 179             if (var[1]==var[0] && var[3]==var[2])
 180               jetSpec->GetCorrelation()->Fill(var,1);
 181           }
 182           else {
 183             // cut at  sigma
 184             if (TMath::Abs(var[1]-var[0]) < sigmax || TMath::Abs(var[3]-var[2]) < sigmay)
 185               jetSpec->GetCorrelation()->Fill(var,func->Eval(var[1],var[3]));
 186           }
 187         }
 188     }
 189
 190
 191   TFile* file = TFile::Open(fileName,"RECREATE");
 192   jetSpec->SaveHistograms();
 193   file->Close();
 194 }
 195
 196 //_____________________________________________________________________________
 197 void StatisticalUncertainties(const char* fileNameGen = "gen_pwg4spec.root", const char* folder = "unfolding", const char* fileNameRec = "rec_pwg4spec.root", const char* fileNameOut = "Uncertainties2DBayesUnf.root")
 198 {
 199   // This function gives the statistical uncertainties due to the 2D Bayesian Unfolding
 200
 201   gSystem->Load("libANALYSIS");
 202   gSystem->Load("libANALYSISalice");
 203   gSystem->Load("libJETAN");
 204   gSystem->Load("libPWG4JetTasks");
 205
 206   AliJetSpectrumUnfolding* jetSpec = new AliJetSpectrumUnfolding(folder, folder);
 207
 208   TFile::Open(fileNameRec);
 209   jetSpec->LoadHistograms();
 210
 211   TFile::Open(fileNameGen);
 212   TH2F* hist = (TH2F*) gFile->Get("unfolding/fGenSpectrum");
 213   jetSpec->SetGenSpectrum(hist);
 214
 215   // create sigma histogram
 216   TH2F* sigma = (TH2F*)(jetSpec->GetGenSpectrum())->Clone("sigma");
 217   sigma->Reset();
 218
 219   THnSparseF* hcorr = (THnSparseF*)(jetSpec->GetCorrelation())->Clone();
 220   TH2F*       htrue = (TH2F*)(jetSpec->GetGenSpectrum())->Clone();
 221   TH2F*       hmeas = (TH2F*)(jetSpec->GetRecSpectrum())->Clone();
 222   TH2F*       hunfo = (TH2F*)(jetSpec->GetUnfSpectrum())->Clone();
 223
 224   Int_t nIterations = 1000;
 225   Float_t binContent;
 226   for(Int_t i=0; i<nIterations; i++)
 227   {
 228     printf("iteration = %d\n", i);
 229     // reset histograms
 230     jetSpec->GetRecSpectrum()->Reset();
 231     jetSpec->GetGenSpectrum()->Reset();
 232     jetSpec->GetCorrelation()->Reset();
 233     jetSpec->GetUnfSpectrum()->Reset();
 234
 235     THnSparseF* tmpcorr = (THnSparseF*)hcorr->Clone("tmpcorr");
 236     TH2F*       tmptrue = (TH2F*)htrue->Clone("tmptrue");
 237
 238     jetSpec->SetGenSpectrum(tmptrue);
 239     jetSpec->SetCorrelation(tmpcorr);
 240
 241     // randomize reconstructed distribution
 242     for (Int_t me=1; me<=hmeas->GetNbinsX(); me++)
 243       for (Int_t mz=1; mz<=hmeas->GetNbinsY(); mz++)
 244       {
 245         binContent = hmeas->GetBinContent(me,mz);
 246         if (binContent>0)
 247         {
 248           TF1* poisson = new TF1("poisson", "TMath::Poisson(x,[0])",binContent*0.25, binContent*1.25);
 249           poisson->SetParameters(binContent,0.);
 250           binContent = poisson->GetRandom();
 251           delete poisson;
 252         }
 253         jetSpec->GetRecSpectrum()->SetBinContent(me,mz, binContent);
 254       }
 255
 256     // unfold
 257     jetSpec->ApplyBayesianMethod(0.2, 20, 0, 0);
 258
 259     // calculate sigma^2
 260     for (Int_t te=1; te<=sigma->GetNbinsX(); te++)
 261       for (Int_t tz=1; tz<=sigma->GetNbinsY(); tz++)
 262       {
 263         if (htrue->GetBinContent(te,tz)!=0)
 264         {
 265           binContent = (jetSpec->GetUnfSpectrum()->GetBinContent(te,tz) -
 266                         htrue->GetBinContent(te,tz))/htrue->GetBinContent(te,tz);
 267           binContent *= binContent;
 268           sigma->SetBinContent(te,tz, binContent + sigma->GetBinContent(te,tz));
 269         }
 270       }
 271   }
 272
 273   // calculate sigma
 274   for (Int_t te=1; te<=sigma->GetNbinsX(); te++)
 275     for (Int_t tz=1; tz<=sigma->GetNbinsY(); tz++)
 276     {
 277       binContent = sigma->GetBinContent(te,tz);
 278       binContent = TMath::Sqrt(binContent/(Float_t)nIterations);
 279       sigma->SetBinContent(te,tz, binContent);
 280     }
 281
 282   const Int_t NRGBs = 5;
 283   const Int_t NCont = 500;
 284
 285   Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
 286   Double_t red[NRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51 };
 287   Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
 288   Double_t blue[NRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00 };
 289   TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
 290   gStyle->SetNumberContours(NCont);
 291   gStyle->SetPalette(1);
 292
 293   new TCanvas();
 294   gPad->SetLogz();
 295   sigma->SetTitle("#sigma((U_{R} - U)/U)");
 296   sigma->SetXTitle("E^{jet} [GeV]");
 297   sigma->SetYTitle("z^{lp}");
 298   sigma->DrawCopy();
 299
 300   TFile* file = TFile::Open(fileNameOut,"RECREATE");
 301   sigma->Write();
 302   file->Close();
 303 }
 304
 305 //_______________________________________________________________________________________________________________
 306 void FillSpecFromFiles(const char* fileNameReal = "histos_pwg4spec.root",const char* fileNameSim = "histos_pwg4spec.root")
 307 {
 308   // This functions avoids problems when the number of bins or the bin limits
 309   // in the histograms of the AliJetSpectrumUnfolding and AliAnalysisTaskJetSpectrum classes
 310   // are different.
 311
 312   gSystem->Load("libANALYSIS");
 313   gSystem->Load("libANALYSISalice");
 314   gSystem->Load("libJETAN");
 315   gSystem->Load("libPWG4JetTasks");
 316
 317   file = new TFile(fileNameSim,"read");
 318   tlist = dynamic_cast<TList*> (file->Get("pwg4spec"));
 319
 320   THnSparseF *fhCorrelation  = 0;
 321   for(int ic = 0;ic<3;++ic){
 322     THnSparseF *hTmp = (THnSparseF*)(tlist->FindObject(Form("fhnCorrelation_%d",ic)));
 323     if(fhCorrelation==0)fhCorrelation =  (THnSparseF*)hTmp->Clone("fhnCorrelation");
 324     else fhCorrelation->Add(hTmp);
 325   }
 326   TH2F *fhEGenZGen = (TH2F*)(tlist->FindObject("fh2EGenZGen"));
 327
 328   AliJetSpectrumUnfolding *jetSpec = new AliJetSpectrumUnfolding("unfolding","unfolding");
 329
 330   // generated jets (true distribution)
 331   for (Int_t te=1; te<=fhEGenZGen->GetNbinsX(); te++)
 332     for (Int_t tz=1; tz<=fhEGenZGen->GetNbinsY(); tz++)
 333     {
 334        Float_t ej = fhEGenZGen->GetXaxis()->GetBinCenter(te);
 335        Float_t  z = fhEGenZGen->GetYaxis()->GetBinCenter(tz);
 336        Int_t bine = jetSpec->GetGenSpectrum()->GetXaxis()->FindBin(ej);
 337        Int_t binz = jetSpec->GetGenSpectrum()->GetYaxis()->FindBin(z);
 338        Float_t cont = jetSpec->GetGenSpectrum()->GetBinContent(bine,binz);
 339        Float_t err  = jetSpec->GetGenSpectrum()->GetBinError(bine,binz);
 340        jetSpec->GetGenSpectrum()->SetBinContent(bine, binz, cont + fhEGenZGen->GetBinContent(te, tz));
 341        jetSpec->GetGenSpectrum()->SetBinError(bine, binz, err + fhEGenZGen->GetBinError(te, tz));
 342     }
 343
 344
 345   Printf("Bins %.0E",jetSpec->GetCorrelation()->GetNbins());
 346   for (Int_t idx=1; idx<=fhCorrelation->GetNbins(); idx++)
 347   {
 348     //printf("idx: %d\n",idx);
 349     Double_t var[4];
 350     Int_t bin[4];
 351     Float_t BinContent = fhCorrelation->GetBinContent(idx, bin);
 352     var[0] = fhCorrelation->GetAxis(0)->GetBinCenter(bin[0]);
 353     var[1] = fhCorrelation->GetAxis(1)->GetBinCenter(bin[1]);
 354     var[2] = fhCorrelation->GetAxis(2)->GetBinCenter(bin[2]);
 355     var[3] = fhCorrelation->GetAxis(3)->GetBinCenter(bin[3]);
 356     bin[0] = jetSpec->GetCorrelation()->GetAxis(0)->FindBin(var[0]);
 357     bin[1] = jetSpec->GetCorrelation()->GetAxis(1)->FindBin(var[1]);
 358     bin[2] = jetSpec->GetCorrelation()->GetAxis(2)->FindBin(var[2]);
 359     bin[3] = jetSpec->GetCorrelation()->GetAxis(3)->FindBin(var[3]);
 360     Float_t cont = jetSpec->GetCorrelation()->GetBinContent(bin);
 361     Float_t err  = jetSpec->GetCorrelation()->GetBinError(bin);
 362     jetSpec->GetCorrelation()->SetBinContent(bin, cont + fhCorrelation->GetBinContent(idx));
 363     jetSpec->GetCorrelation()->SetBinError(bin, err + fhCorrelation->GetBinError(idx));
 364   }
 365   // need number of entries for correct drawing
 366   jetSpec->GetCorrelation()->SetEntries(fhCorrelation->GetEntries());
 367
 368
 369   file = TFile::Open("gen_pwg4spec.root", "RECREATE");
 370   jetSpec->SaveHistograms();
 371   Printf("Bins %.0E",jetSpec->GetCorrelation()->GetNbins());
 372   file->Close();
 373   jetSpec->GetGenSpectrum()->Reset();
 374   printf("true distribution has been set\n");
 375
 376   // reconstructed jets (measured distribution)
 377
 378
 379   // Response function
 380   jetSpec->GetCorrelation()->Reset();
 381   jetSpec->GetCorrelation()->Sumw2();
 382   jetSpec->GetGenSpectrum()->Reset();
 383   jetSpec->GetRecSpectrum()->Reset();
 384
 385
 386   file = new TFile(fileNameReal,"read");
 387   tlist = dynamic_cast<TList*> (file->Get("pwg4spec"));
 388
 389
 390   TH2F *fhERecZRec = (TH2F*)(tlist->FindObject("fh2ERecZRec"));
 391
 392   for (Int_t me=1; me<=fhERecZRec->GetNbinsX(); me++)
 393     for (Int_t mz=1; mz<=fhERecZRec->GetNbinsY(); mz++)
 394     {
 395        Float_t erec = fhERecZRec->GetXaxis()->GetBinCenter(me);
 396        Float_t   zm = fhERecZRec->GetYaxis()->GetBinCenter(mz);
 397        Int_t bine   = jetSpec->GetRecSpectrum()->GetXaxis()->FindBin(erec);
 398        Int_t binz   = jetSpec->GetRecSpectrum()->GetYaxis()->FindBin(zm);
 399        Float_t cont = jetSpec->GetRecSpectrum()->GetBinContent(bine, binz);
 400        Float_t err  = jetSpec->GetRecSpectrum()->GetBinError(bine, binz);
 401        jetSpec->GetRecSpectrum()->SetBinContent(bine, binz, cont + fhERecZRec->GetBinContent(me, mz));
 402        jetSpec->GetRecSpectrum()->SetBinError(bine, binz, err + fhERecZRec->GetBinError(me, mz));
 403     }
 404
 405
 406   // for control again, but now from the rec file
 407   // generated jets (true distribution)
 408   TH2F *fhEGenZGen = (TH2F*)(tlist->FindObject("fh2EGenZGen"));
 409   for (Int_t te=1; te<=fhEGenZGen->GetNbinsX(); te++)
 410     for (Int_t tz=1; tz<=fhEGenZGen->GetNbinsY(); tz++)
 411     {
 412        Float_t ej = fhEGenZGen->GetXaxis()->GetBinCenter(te);
 413        Float_t  z = fhEGenZGen->GetYaxis()->GetBinCenter(tz);
 414        Int_t bine = jetSpec->GetGenSpectrum()->GetXaxis()->FindBin(ej);
 415        Int_t binz = jetSpec->GetGenSpectrum()->GetYaxis()->FindBin(z);
 416        Float_t cont = jetSpec->GetGenSpectrum()->GetBinContent(bine,binz);
 417        Float_t err  = jetSpec->GetGenSpectrum()->GetBinError(bine,binz);
 418        jetSpec->GetGenSpectrum()->SetBinContent(bine, binz, cont + fhEGenZGen->GetBinContent(te, tz));
 419        jetSpec->GetGenSpectrum()->SetBinError(bine, binz, err + fhEGenZGen->GetBinError(te, tz));
 420     }
 421
 422
 423   file = TFile::Open("rec_pwg4spec.root", "RECREATE");
 424   jetSpec->SaveHistograms();
 425   file->Close();
 426
 427   printf("reconstructed distribution has been set\n");
 428   printf("response map has been set\n");
 429
 430 }
 431
 432 void correct(){
 433   // simple steering to correct a given distribution;
 434   loadlibs();
 435   // rec and gen
 436   //  FillSpecFromFiles("pwg4spec_15-50_all.root","pwg4spec_allpt.root");
 437   FillSpecFromFiles("pwg4spec_allpt.root","pwg4spec_allpt.root");
 438
 439   char name[100];
 440   sprintf(name, "unfolded_pwg4spec.root");
 441
 442   unfold("gen_pwg4spec.root", "unfolding", "rec_pwg4spec.root", name);
 443   //draw(name, "unfolding", 1);
 444
 445 }
 446
 447 void mergeJetAnaOutput(){
 448   // This is used to merge the analysis-output from different
 449   // data samples/pt_hard bins
 450   // in case the eventweigth was set to xsection/ntrials already, this
 451   // is not needed. Both methods only work in case we do not mix different
 452   // pt_hard bins, and do not have overlapping bins
 453
 454   const Int_t nBins = 2;
 455   // LHC08q jetjet100: Mean = 1.42483e-03, RMS = 6.642e-05
 456   // LHC08r jetjet50: Mean = 2.44068e-02, RMS = 1.144e-03
 457   // LHC08v jetjet15-50: Mean = 2.168291 , RMS = 7.119e-02
 458   // const Float_t xsection[nBins] = {2.168291,2.44068e-02};
 459   Float_t xsection[nBins] = {0,0};
 460   Float_t nTrials[nBins] = {0,0};
 461   Float_t sf[nBins] = {0,0};
 462
 463   const char *cFile[nBins] = {"pwg4spec_0000000-1000000_LHC08v_jetjet15-50.root",
 464                               "pwg4spec_0000000-1000000_LHC08r_jetjet50.root"};
 465
 466
 467   TList *lIn[nBins];
 468   TFile *fIn[nBins];
 469   for(int ib = 0;ib < nBins;++ib){
 470     fIn[ib] = TFile::Open(cFile[ib]);
 471     lIn[ib] = (TList*)fIn[ib]->Get("pwg4spec");
 472     TH1* hTrials = (TH1F*)lIn[ib]->FindObject("fh1PtHard_Trials");
 473     TProfile* hXsec = (TProfile*)lIn[ib]->FindObject("fh1Xsec");
 474     xsection[ib] = hXsec->GetBinContent(1);
 475     nTrials[ib] = hTrials->Integral();
 476     sf[ib] = xsection[ib]/ nTrials[ib];
 477   }
 478
 479   TFile *fOut = new TFile("pwg4spec_allpt.root","RECREATE");
 480   TList *lOut = new TList();
 481   lOut->SetName(lIn[0]->GetName());
 482   // for the start scale all...
 483   for(int ie = 0; ie < lIn[0]->GetEntries();++ie){
 484     TH1 *h1Add = 0;
 485     THnSparse *hnAdd = 0;
 486     Printf("%d: %s",ie, lIn[0]->At(ie)->GetName());
 487     for(int ib = 0;ib < nBins;++ib){
 488       // dynamic cast does not work with cint
 489       TObject *h = lIn[ib]->At(ie);
 490       if(h->InheritsFrom("TH1")){
 491         Printf("ib %d",ib);
 492         TH1 *h1 = (TH1*)h;
 493         if(ib==0){
 494           h1Add = (TH1*)h1->Clone(h1->GetName());
 495           h1Add->Scale(sf[ib]);
 496         }
 497         else{
 498           h1Add->Add(h1,sf[ib]);
 499         }
 500       }
 501       else if(h->InheritsFrom("THnSparse")){
 502         Printf("ib %d",ib);
 503         THnSparse *hn = (THnSparse*)h;
 504         if(ib==0){
 505           hnAdd = (THnSparse*)hn->Clone(hn->GetName());
 506           hnAdd->Scale(sf[ib]);
 507         }
 508         else{
 509           hnAdd->Add(hn,sf[ib]);
 510         }
 511       }
 512
 513
 514     }// ib
 515     if(h1Add)lOut->Add(h1Add);
 516     else if(hnAdd)lOut->Add(hnAdd);
 517   }
 518   fOut->cd();
 519   lOut->Write(lOut->GetName(),TObject::kSingleKey);
 520   fOut->Close();
 521
 522
 523
 524 }