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        // merging of bins happens here!
 341        jetSpec->GetGenSpectrum()->SetBinContent(bine, binz, cont + fhEGenZGen->GetBinContent(te, tz));
 342        jetSpec->GetGenSpectrum()->SetBinError(bine, binz, err + fhEGenZGen->GetBinError(te, tz));
 343     }
 344
 345
 346   Printf("Bins %.0E",jetSpec->GetCorrelation()->GetNbins());
 347   for (Int_t idx=1; idx<=fhCorrelation->GetNbins(); idx++)
 348   {
 349     //printf("idx: %d\n",idx);
 350     Double_t var[4];
 351     Int_t bin[4];
 352     Float_t BinContent = fhCorrelation->GetBinContent(idx, bin);
 353     var[0] = fhCorrelation->GetAxis(0)->GetBinCenter(bin[0]);
 354     var[1] = fhCorrelation->GetAxis(1)->GetBinCenter(bin[1]);
 355     var[2] = fhCorrelation->GetAxis(2)->GetBinCenter(bin[2]);
 356     var[3] = fhCorrelation->GetAxis(3)->GetBinCenter(bin[3]);
 357     bin[0] = jetSpec->GetCorrelation()->GetAxis(0)->FindBin(var[0]);
 358     bin[1] = jetSpec->GetCorrelation()->GetAxis(1)->FindBin(var[1]);
 359     bin[2] = jetSpec->GetCorrelation()->GetAxis(2)->FindBin(var[2]);
 360     bin[3] = jetSpec->GetCorrelation()->GetAxis(3)->FindBin(var[3]);
 361     Float_t cont = jetSpec->GetCorrelation()->GetBinContent(bin);
 362     Float_t err  = jetSpec->GetCorrelation()->GetBinError(bin);
 363     // merging of bins happens here!
 364     jetSpec->GetCorrelation()->SetBinContent(bin, cont + fhCorrelation->GetBinContent(idx));
 365     jetSpec->GetCorrelation()->SetBinError(bin, err + fhCorrelation->GetBinError(idx));
 366   }
 367   // need number of entries for correct drawing
 368   jetSpec->GetCorrelation()->SetEntries(fhCorrelation->GetEntries());
 369
 370
 371   file = TFile::Open("gen_pwg4spec.root", "RECREATE");
 372   jetSpec->SaveHistograms();
 373   Printf("Bins %.0E",jetSpec->GetCorrelation()->GetNbins());
 374   file->Close();
 375   jetSpec->GetGenSpectrum()->Reset();
 376   printf("true distribution has been set\n");
 377
 378   // reconstructed jets (measured distribution)
 379
 380
 381   // Response function
 382   jetSpec->GetCorrelation()->Reset();
 383   jetSpec->GetCorrelation()->Sumw2();
 384   jetSpec->GetGenSpectrum()->Reset();
 385   jetSpec->GetRecSpectrum()->Reset();
 386
 387
 388   file = new TFile(fileNameReal,"read");
 389   tlist = dynamic_cast<TList*> (file->Get("pwg4spec"));
 390
 391
 392   TH2F *fhERecZRec = (TH2F*)(tlist->FindObject("fh2ERecZRec"));
 393
 394   for (Int_t me=1; me<=fhERecZRec->GetNbinsX(); me++)
 395     for (Int_t mz=1; mz<=fhERecZRec->GetNbinsY(); mz++)
 396     {
 397        Float_t erec = fhERecZRec->GetXaxis()->GetBinCenter(me);
 398        Float_t   zm = fhERecZRec->GetYaxis()->GetBinCenter(mz);
 399        Int_t bine   = jetSpec->GetRecSpectrum()->GetXaxis()->FindBin(erec);
 400        Int_t binz   = jetSpec->GetRecSpectrum()->GetYaxis()->FindBin(zm);
 401        Float_t cont = jetSpec->GetRecSpectrum()->GetBinContent(bine, binz);
 402        Float_t err  = jetSpec->GetRecSpectrum()->GetBinError(bine, binz);
 403        jetSpec->GetRecSpectrum()->SetBinContent(bine, binz, cont + fhERecZRec->GetBinContent(me, mz));
 404        jetSpec->GetRecSpectrum()->SetBinError(bine, binz, err + fhERecZRec->GetBinError(me, mz));
 405     }
 406
 407
 408   // for control again, but now from the rec file
 409   // generated jets (true distribution)
 410   TH2F *fhEGenZGen = (TH2F*)(tlist->FindObject("fh2EGenZGen"));
 411   for (Int_t te=1; te<=fhEGenZGen->GetNbinsX(); te++)
 412     for (Int_t tz=1; tz<=fhEGenZGen->GetNbinsY(); tz++)
 413     {
 414        Float_t ej = fhEGenZGen->GetXaxis()->GetBinCenter(te);
 415        Float_t  z = fhEGenZGen->GetYaxis()->GetBinCenter(tz);
 416        Int_t bine = jetSpec->GetGenSpectrum()->GetXaxis()->FindBin(ej);
 417        Int_t binz = jetSpec->GetGenSpectrum()->GetYaxis()->FindBin(z);
 418        Float_t cont = jetSpec->GetGenSpectrum()->GetBinContent(bine,binz);
 419        Float_t err  = jetSpec->GetGenSpectrum()->GetBinError(bine,binz);
 420        jetSpec->GetGenSpectrum()->SetBinContent(bine, binz, cont + fhEGenZGen->GetBinContent(te, tz));
 421        jetSpec->GetGenSpectrum()->SetBinError(bine, binz, err + fhEGenZGen->GetBinError(te, tz));
 422     }
 423
 424
 425   file = TFile::Open("rec_pwg4spec.root", "RECREATE");
 426   jetSpec->SaveHistograms();
 427   file->Close();
 428
 429   printf("reconstructed distribution has been set\n");
 430   printf("response map has been set\n");
 431
 432 }
 433
 434 void correct(){
 435   // simple steering to correct a given distribution;
 436   loadlibs();
 437   // rec and gen
 438   //  FillSpecFromFiles("pwg4spec_15-50_all.root","pwg4spec_allpt.root");
 439   FillSpecFromFiles("pwg4spec_allpt_16.root","pwg4spec_allpt_16.root");
 440
 441   char name[100];
 442   sprintf(name, "unfolded_pwg4spec_16.root");
 443
 444   unfold("gen_pwg4spec.root", "unfolding", "rec_pwg4spec.root", name);
 445   //draw(name, "unfolding", 1);
 446
 447 }
 448
 449 void mergeJetAnaOutput(){
 450   // This is used to merge the analysis-output from different
 451   // data samples/pt_hard bins
 452   // in case the eventweigth was set to xsection/ntrials already, this
 453   // is not needed. Both methods only work in case we do not mix different
 454   // pt_hard bins, and do not have overlapping bins
 455
 456   const Int_t nBins = 2;
 457   // LHC08q jetjet100: Mean = 1.42483e-03, RMS = 6.642e-05
 458   // LHC08r jetjet50: Mean = 2.44068e-02, RMS = 1.144e-03
 459   // LHC08v jetjet15-50: Mean = 2.168291 , RMS = 7.119e-02
 460   // const Float_t xsection[nBins] = {2.168291,2.44068e-02};
 461   Float_t xsection[nBins] = {0,0};
 462   Float_t nTrials[nBins] = {0,0};
 463   Float_t sf[nBins] = {0,0};
 464
 465   const char *cFile[nBins] = {"pwg4spec_0000000-1000000_LHC08v_jetjet15-50.root",
 466                               "pwg4spec_0000000-1000000_LHC08r_jetjet50.root"};
 467
 468
 469   TList *lIn[nBins];
 470   TFile *fIn[nBins];
 471   for(int ib = 0;ib < nBins;++ib){
 472     fIn[ib] = TFile::Open(cFile[ib]);
 473     lIn[ib] = (TList*)fIn[ib]->Get("pwg4spec");
 474     TH1* hTrials = (TH1F*)lIn[ib]->FindObject("fh1PtHard_Trials");
 475     TProfile* hXsec = (TProfile*)lIn[ib]->FindObject("fh1Xsec");
 476     xsection[ib] = hXsec->GetBinContent(1);
 477     nTrials[ib] = hTrials->Integral();
 478     sf[ib] = xsection[ib]/ nTrials[ib];
 479   }
 480
 481   TFile *fOut = new TFile("pwg4spec_allpt.root","RECREATE");
 482   TList *lOut = new TList();
 483   lOut->SetName(lIn[0]->GetName());
 484   // for the start scale all...
 485   for(int ie = 0; ie < lIn[0]->GetEntries();++ie){
 486     TH1 *h1Add = 0;
 487     THnSparse *hnAdd = 0;
 488     Printf("%d: %s",ie, lIn[0]->At(ie)->GetName());
 489     for(int ib = 0;ib < nBins;++ib){
 490       // dynamic cast does not work with cint
 491       TObject *h = lIn[ib]->At(ie);
 492       if(h->InheritsFrom("TH1")){
 493         Printf("ib %d",ib);
 494         TH1 *h1 = (TH1*)h;
 495         if(ib==0){
 496           h1Add = (TH1*)h1->Clone(h1->GetName());
 497           h1Add->Scale(sf[ib]);
 498         }
 499         else{
 500           h1Add->Add(h1,sf[ib]);
 501         }
 502       }
 503       else if(h->InheritsFrom("THnSparse")){
 504         Printf("ib %d",ib);
 505         THnSparse *hn = (THnSparse*)h;
 506         if(ib==0){
 507           hnAdd = (THnSparse*)hn->Clone(hn->GetName());
 508           hnAdd->Scale(sf[ib]);
 509         }
 510         else{
 511           hnAdd->Add(hn,sf[ib]);
 512         }
 513       }
 514
 515
 516     }// ib
 517     if(h1Add)lOut->Add(h1Add);
 518     else if(hnAdd)lOut->Add(hnAdd);
 519   }
 520   fOut->cd();
 521   lOut->Write(lOut->GetName(),TObject::kSingleKey);
 522   fOut->Close();
 523
 524
 525
 526 }