]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWG4/macros/runJetSpectrumUnfolding.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Cleaning up files, added documentation. Added ntuple for v0s.
[u/mrichter/AliRoot.git] / PWG4 / macros / runJetSpectrumUnfolding.C
CommitLineData
6d75bdb8 1//
2// script with functions to use AliJetSpectrumUnfolding class
3//
4
5
6void 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
26void 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//________________________________________________________________________________________________________________
74void 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
df65bddb 83 TFile::Open(fileNameGen);
6d75bdb8 84 jetSpec->LoadHistograms();
85
df65bddb 86 TFile::Open(fileNameRec);
87 TH2F* hist = (TH2F*) gFile->Get("unfolding/fRecSpectrum");
88 jetSpec->SetRecSpectrum(hist);
4dbfdecc 89 hist = (TH2F*) gFile->Get("unfolding/fGenSpectrum");
90 if(hist->GetEntries()>0)jetSpec->SetGenSpectrum(hist);
6d75bdb8 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//___________________________________________________________________________
101void 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//_____________________________________________________________________________________________
150void 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//_____________________________________________________________________________
197void 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//_______________________________________________________________________________________________________________
df65bddb 306void FillSpecFromFiles(const char* fileNameReal = "histos_pwg4spec.root",const char* fileNameSim = "histos_pwg4spec.root")
6d75bdb8 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
df65bddb 317 file = new TFile(fileNameSim,"read");
6d75bdb8 318 tlist = dynamic_cast<TList*> (file->Get("pwg4spec"));
6d75bdb8 319
df65bddb 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"));
6d75bdb8 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 }
6d75bdb8 343
6d75bdb8 344
df65bddb 345 Printf("Bins %.0E",jetSpec->GetCorrelation()->GetNbins());
6d75bdb8 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 }
df65bddb 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 }
6d75bdb8 404
4dbfdecc 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
6d75bdb8 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
6d75bdb8 432void correct(){
433 // simple steering to correct a given distribution;
434 loadlibs();
4dbfdecc 435 // rec and gen
436 // FillSpecFromFiles("pwg4spec_15-50_all.root","pwg4spec_allpt.root");
437 FillSpecFromFiles("pwg4spec_allpt.root","pwg4spec_allpt.root");
6d75bdb8 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}
4dbfdecc 446
447void 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 cse 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 nTrials[nBins] = {0,0};
460 Float_t sf[nBins] = {0,0};
461
462 const char *cFile[nBins] = {"pwg4spec_15-50_all.root","pwg4spec_50_all.root"};
463
464
465 TList *lIn[nBins];
466 TFile *fIn[nBins];
467 for(int ib = 0;ib < nBins;++ib){
468 fIn[ib] = TFile::Open(cFile[ib]);
469 lIn[ib] = (TList*)fIn[ib]->Get("pwg4spec");
470 TH1* hTrials = (TH1F*)lIn[ib]->FindObject("fh1PtHard_Trials");
471 nTrials[ib] = hTrials->Integral();
472 sf[ib] = xsection[ib]/ nTrials[ib];
473 }
474
475 TFile *fOut = new TFile("pwg4spec_allpt.root","RECREATE");
476 TList *lOut = new TList();
477 lOut->SetName(lIn[0]->GetName());
478 // for the start scale all...
479 for(int ie = 0; ie < lIn[0]->GetEntries();++ie){
480 TH1 *h1Add = 0;
481 THnSparse *hnAdd = 0;
482 Printf("%d: %s",ie, lIn[0]->At(ie)->GetName());
483 for(int ib = 0;ib < nBins;++ib){
484 // dynamic cast does not work with cint
485 TObject *h = lIn[ib]->At(ie);
486 if(h->InheritsFrom("TH1")){
487 Printf("ib %d",ib);
488 TH1 *h1 = (TH1*)h;
489 if(ib==0){
490 h1Add = (TH1*)h1->Clone(h1->GetName());
491 h1Add->Scale(sf[ib]);
492 }
493 else{
494 h1Add->Add(h1,sf[ib]);
495 }
496 }
497 else if(h->InheritsFrom("THnSparse")){
498 Printf("ib %d",ib);
499 THnSparse *hn = (THnSparse*)h;
500 if(ib==0){
501 hnAdd = (THnSparse*)hn->Clone(hn->GetName());
502 hnAdd->Scale(sf[ib]);
503 }
504 else{
505 hnAdd->Add(hn,sf[ib]);
506 }
507 }
508
509
510 }// ib
511 if(h1Add)lOut->Add(h1Add);
512 else if(hnAdd)lOut->Add(hnAdd);
513 }
514 fOut->cd();
515 lOut->Write(lOut->GetName(),TObject::kSingleKey);
516 fOut->Close();
517
518
519
520}