PWGHF/vertexingHF/macros/CombineFeedDownMCSubtractionMethodsUncertainties.C

   1 #if !defined(__CINT__) || defined(__MAKECINT__)
   2 #include "TFile.h"
   3 #include "TH1.h"
   4 #include "TH1D.h"
   5 #include "TH2.h"
   6 #include "TH2F.h"
   7 #include "TGraphAsymmErrors.h"
   8 #include "TCanvas.h"
   9 #include "TLegend.h"
  10 #include "TMath.h"
  11 #include "TROOT.h"
  12 #include "TStyle.h"
  13 #include "AliHFSystErr.h"
  14 #include <Riostream.h>
  15 #endif
  16
  17 /* $Id$ */
  18
  19
  20 //_________________________________________________________________________________________
  21 //
  22 //  Macro to combine the the MonteCarlo B feed-down subtraction uncertainties
  23 //
  24 //   Take as input the output files from the HFPtSpectrum class
  25 //    from both fc & Nb subtraction methods and combine the uncertainties.
  26 //   The final central value is set as the one from the Nb-method.
  27 //   The final uncertainties are defined as the envelope of both fc & Nb
  28 //      uncertainties with respect to the new central-value.
  29 //   The final global uncertainties are also defined and a preliminary drawing done.
  30 //
  31 //
  32 //   Usage parameters:
  33 //      1. HFPtSpectrum fc subtraction file
  34 //      2. HFPtSpectrum Nb subtraction file
  35 //      3. Output file name
  36 //      4. FONLL theoretical predictions file to draw on top
  37 //      5. Decay channel as defined in the AliHFSystErr class
  38 //
  39 //_________________________________________________________________________________________
  40
  41 enum centrality{ kpp7, kpp276, k010, k020, k2040, k4060, k6080, k4080, k80100 };
  42
  43 void CombineFeedDownMCSubtractionMethodsUncertainties(const char *fcfilename="HFPtSpectrum_D0Kpi_method1_221110_newnorm.root",
  44                                                       const char *nbfilename="HFPtSpectrum_D0Kpi_method2_221110_newnorm.root",
  45                                                       const char *outfilename="HFPtSpectrum_D0Kpi_combinedFD.root",
  46                                                       const char *thfilename="D0DplusDstarPredictions_y05.root",
  47                                                       Int_t decay=1, Int_t centrality=kpp7)
  48 {
  49
  50   //
  51   // Get fc file inputs
  52   TFile * fcfile = new TFile(fcfilename,"read");
  53   TH1D * histoSigmaCorrFc = (TH1D*)fcfile->Get("histoSigmaCorr");
  54   histoSigmaCorrFc->SetNameTitle("histoSigmaCorrFc","histoSigmaCorrFc");
  55   TGraphAsymmErrors * gSigmaCorrFc = (TGraphAsymmErrors*)fcfile->Get("gSigmaCorr");
  56   gSigmaCorrFc->SetNameTitle("gSigmaCorrFc","gSigmaCorrFc");
  57   TGraphAsymmErrors * gSigmaCorrConservativeFc = (TGraphAsymmErrors*)fcfile->Get("gSigmaCorrConservative");
  58   gSigmaCorrConservativeFc->SetNameTitle("gSigmaCorrConservativeFc","Cross section (fc prompt fraction)");
  59   TGraphAsymmErrors * gFcConservativeFc = (TGraphAsymmErrors*)fcfile->Get("gFcConservative");
  60   gFcConservativeFc->SetNameTitle("gFcConservativeFc","fc prompt fraction");
  61
  62   //
  63   // Get Nb file inputs
  64   TFile * nbfile = new TFile(nbfilename,"read");
  65   TH1D * histoSigmaCorrNb = (TH1D*)nbfile->Get("histoSigmaCorr");
  66   histoSigmaCorrNb->SetNameTitle("histoSigmaCorrNb","histoSigmaCorrNb");
  67   TGraphAsymmErrors * gSigmaCorrNb = (TGraphAsymmErrors*)nbfile->Get("gSigmaCorr");
  68   gSigmaCorrNb->SetNameTitle("gSigmaCorrNb","gSigmaCorrNb");
  69   TGraphAsymmErrors * gSigmaCorrConservativeNb = (TGraphAsymmErrors*)nbfile->Get("gSigmaCorrConservative");
  70   gSigmaCorrConservativeNb->SetNameTitle("gSigmaCorrConservativeNb","Cross section (Nb prompt fraction)");
  71   TGraphAsymmErrors * gFcConservativeNb = (TGraphAsymmErrors*)nbfile->Get("gFcConservative");
  72   gFcConservativeNb->SetNameTitle("gFcConservativeNb","Nb prompt fraction");
  73
  74   //
  75   // Get the predictions input
  76   TFile *thfile = new TFile(thfilename,"read");
  77   TGraphAsymmErrors * thD0KpifromBprediction = (TGraphAsymmErrors*)thfile->Get("D0Kpiprediction");
  78   TGraphAsymmErrors * thDpluskpipiprediction = (TGraphAsymmErrors*)thfile->Get("Dpluskpipiprediction");
  79   TGraphAsymmErrors * thDstarD0piprediction = (TGraphAsymmErrors*)thfile->Get("DstarD0piprediction");
  80   TGraphAsymmErrors * thDsKKpiprediction = (TGraphAsymmErrors*)thfile->Get("DsKkpiprediction");
  81
  82   thD0KpifromBprediction->SetLineColor(4);
  83   thD0KpifromBprediction->SetFillColor(kAzure+9);
  84   thDpluskpipiprediction->SetLineColor(4);
  85   thDpluskpipiprediction->SetFillColor(kAzure+9);
  86   thDstarD0piprediction->SetLineColor(4);
  87   thDstarD0piprediction->SetFillColor(kAzure+9);
  88   thDsKKpiprediction->SetLineColor(4);
  89   thDsKKpiprediction->SetFillColor(kAzure+9);
  90
  91   //
  92   // Get the spectra bins & limits
  93   Int_t nbins = histoSigmaCorrFc->GetNbinsX();
  94   Double_t *limits = new Double_t[nbins+1];
  95   Double_t xlow=0., binwidth=0.;
  96   for (Int_t i=1; i<=nbins; i++) {
  97     binwidth = histoSigmaCorrFc->GetBinWidth(i);
  98     xlow = histoSigmaCorrFc->GetBinLowEdge(i);
  99     limits[i-1] = xlow;
 100   }
 101   limits[nbins] = xlow + binwidth;
 102
 103
 104   //
 105   // Define a new histogram with the real-data reconstructed spectrum binning
 106   //   they will be filled with central value equal to the Nb result
 107   //   and uncertainties taken from the envelope of the result uncertainties
 108   //   The systematical unc. (but FD) will also be re-calculated
 109   //
 110   TH1D * histoSigmaCorr = new TH1D("histoSigmaCorr","corrected cross-section (combined fc and Nb MC feed-down subtraction)",nbins,limits);
 111   TGraphAsymmErrors * gSigmaCorr = new TGraphAsymmErrors(nbins+1);
 112   gSigmaCorr->SetNameTitle("gSigmaCorr","gSigmaCorr (combined fc and Nb MC FD)");
 113   TGraphAsymmErrors * gFcCorrConservative = new TGraphAsymmErrors(nbins+1);
 114   gFcCorrConservative->SetNameTitle("gFcCorrConservative","Combined prompt fraction");
 115   TGraphAsymmErrors * gSigmaCorrConservative = new TGraphAsymmErrors(nbins+1);
 116   gSigmaCorrConservative->SetNameTitle("gSigmaCorrConservative","Cross section (combined prompt fraction)");
 117   TGraphAsymmErrors * gSigmaCorrConservativePC = new TGraphAsymmErrors(nbins+1);
 118   gSigmaCorrConservativePC->SetNameTitle("gSigmaCorrConservativePC","Conservative gSigmaCorr (combined fc and Nb MC FD) in percentages [for drawing with AliHFSystErr]");
 119
 120   //
 121   // Call the systematics uncertainty class for a given decay
 122   //   will help to compute the systematical unc. (but FD)
 123   AliHFSystErr systematics;
 124   if( centrality==kpp276 ) {
 125     systematics.SetIsLowEnergy(true);
 126   } else if( centrality!=kpp7 )  {
 127     systematics.SetCollisionType(1);
 128     if ( centrality == k020 ) {
 129       systematics.SetCentrality("020");
 130     }
 131     else if ( centrality == k4080 ) {
 132       systematics.SetCentrality("4080");
 133     }
 134     else {
 135       cout << " Systematics not yet implemented " << endl;
 136       return;
 137     }
 138   }
 139   else { systematics.SetCollisionType(0); }
 140   systematics.Init(decay);
 141
 142   //
 143   // Loop on all the bins to do the calculations
 144   //
 145   Double_t pt=0., average = 0., averageStatUnc=0., avErrx=0., avErryl=0., avErryh=0., avErryfdl=0., avErryfdh=0.;
 146   Double_t avErrylPC=0., avErryhPC=0., avErryfdlPC=0., avErryfdhPC=0.;
 147   Double_t valFc = 0., valFcErrstat=0., valFcErrx=0., valFcErryl=0., valFcErryh=0., valFcErryfdl=0., valFcErryfdh=0.;
 148   Double_t valNb = 0., valNbErrstat=0., valNbErrx=0., valNbErryl=0., valNbErryh=0., valNbErryfdl=0., valNbErryfdh=0.;
 149   Double_t corrfd = 0., corrfdl=0., corrfdh=0.;
 150   //
 151   for(Int_t ibin=1; ibin<=nbins; ibin++){
 152
 153     // Get input values from fc method
 154     valFc = histoSigmaCorrFc->GetBinContent(ibin);
 155     pt = histoSigmaCorrFc->GetBinCenter(ibin);
 156     valFcErrstat = histoSigmaCorrFc->GetBinError(ibin);
 157     Double_t value =0., ptt=0.;
 158     gSigmaCorrConservativeFc->GetPoint(ibin,ptt,value);
 159     if (value<=0.) continue;
 160     if ( TMath::Abs(valFc-value)>0.1 || TMath::Abs(pt-ptt)>0.1 )
 161       cout << "Hey you ! There might be a problem with the fc input file, please, have a look !" << endl;
 162     valFcErrx = gSigmaCorrFc->GetErrorXlow(ibin);
 163     valFcErryl = gSigmaCorrFc->GetErrorYlow(ibin);
 164     valFcErryh = gSigmaCorrFc->GetErrorYhigh(ibin);
 165     valFcErryfdl = TMath::Abs( gSigmaCorrConservativeFc->GetErrorYlow(ibin) );
 166     valFcErryfdh = TMath::Abs( gSigmaCorrConservativeFc->GetErrorYhigh(ibin) );
 167     Double_t valfdFc = 0., x=0.;
 168     gFcConservativeFc->GetPoint(ibin,x,valfdFc);
 169     Double_t valfdFch = gFcConservativeFc->GetErrorYhigh(ibin);
 170     Double_t valfdFcl = gFcConservativeFc->GetErrorYlow(ibin);
 171
 172     // Get input values from Nb method
 173     valNb = histoSigmaCorrNb->GetBinContent(ibin);
 174     pt = histoSigmaCorrNb->GetBinCenter(ibin);
 175     valNbErrstat = histoSigmaCorrNb->GetBinError(ibin);
 176     gSigmaCorrConservativeNb->GetPoint(ibin,ptt,value);
 177     if ( TMath::Abs(valNb-value)>0.1 || TMath::Abs(pt-ptt)>0.1 )
 178       cout << "Hey you ! There might be a problem with the Nb input file, please, have a look !" << endl;
 179     valNbErrx = gSigmaCorrNb->GetErrorXlow(ibin);
 180     valNbErryl = gSigmaCorrNb->GetErrorYlow(ibin);
 181     valNbErryh = gSigmaCorrNb->GetErrorYhigh(ibin);
 182     valNbErryfdl = gSigmaCorrConservativeNb->GetErrorYlow(ibin);
 183     valNbErryfdh = gSigmaCorrConservativeNb->GetErrorYhigh(ibin);
 184     Double_t valfdNb = 0.;
 185     gFcConservativeNb->GetPoint(ibin,x,valfdNb);
 186     Double_t valfdNbh = gFcConservativeNb->GetErrorYhigh(ibin);
 187     Double_t valfdNbl = gFcConservativeNb->GetErrorYlow(ibin);
 188
 189
 190     // Compute the FD combined value
 191     //    average = valNb
 192     average = valNb ;
 193     corrfd = valfdNb;
 194     avErrx = valFcErrx;
 195     if ( TMath::Abs( valFcErrx - valNbErrx ) > 0.1 )
 196       cout << "Hey you ! There might be consistency problem with the fc & Nb input files, please, have a look !" << endl;
 197     averageStatUnc = valNbErrstat ;
 198 //     cout << " pt=" << pt << ", average="<<average<<endl;
 199 //     cout << "   stat unc (pc)=" << averageStatUnc/average << ", stat-fc (pc)="<<(valFcErrstat/valFc) << ", stat-Nb (pc)="<<(valNbErrstat/valNb)<<endl;
 200
 201     // now estimate the new feed-down combined uncertainties
 202     Double_t minimum[2] = { (valFc - valFcErryfdl), (valNb - valNbErryfdl) };
 203     Double_t maximum[2] = { (valFc + valFcErryfdh), (valNb + valNbErryfdh) };
 204     avErryfdl = average - TMath::MinElement(2,minimum);
 205     avErryfdh = TMath::MaxElement(2,maximum) - average;
 206     avErryfdlPC = avErryfdl / average ; // in percentage
 207     avErryfdhPC = avErryfdh / average ; // in percentage
 208 //     cout << " fc : val " << valFc << " + " << valFcErryfdh <<" - " << valFcErryfdl <<endl;
 209 //     cout << " Nb : val " << valNb << " + " << valNbErryfdh <<" - " << valNbErryfdl <<endl;
 210 //     cout << " fc  & Nb: val " << average << " + " << avErryfdh <<" - " << avErryfdl <<endl;
 211     Double_t minimumfc[2] = { (valfdNb - valfdNbl), (valfdFc - valfdFcl) };
 212     Double_t maximumfc[2] = { (valfdNb + valfdNbh), (valfdFc + valfdFch) };
 213     corrfdl = corrfd - TMath::MinElement(2,minimumfc);
 214     corrfdh = TMath::MaxElement(2,maximumfc) - corrfd;
 215
 216
 217     // compute the global systematics
 218     avErrylPC = systematics.GetTotalSystErr(pt,avErryfdlPC); // in percentage
 219     avErryhPC = systematics.GetTotalSystErr(pt,avErryfdhPC); // in percentage
 220     avErryl = avErrylPC * average ;
 221     avErryh = avErryhPC * average ;
 222 //     cout << "   syst av-l="<<avErryl<<", av-h="<<avErryh<<endl;
 223 //     cout << "     fd-l-pc="<<avErryfdlPC<<", fd-h-pc="<<avErryfdhPC<<", syst err(no fd)-pc="<<systematics.GetTotalSystErr(pt)<<", av-l-pc="<<avErrylPC<<", av-h-pc="<<avErryhPC<<endl;
 224
 225     // fill in the histos and TGraphs
 226     //   fill them only when for non empty bins
 227     if ( average > 0.1 ) {
 228       histoSigmaCorr->SetBinContent(ibin,average);
 229       histoSigmaCorr->SetBinError(ibin,averageStatUnc);
 230       gSigmaCorr->SetPoint(ibin,pt,average);
 231       gSigmaCorr->SetPointError(ibin,valFcErrx,valFcErrx,avErryl,avErryh);
 232       gSigmaCorrConservative->SetPoint(ibin,pt,average);
 233       gSigmaCorrConservative->SetPointError(ibin,valFcErrx,valFcErrx,avErryfdl,avErryfdh);
 234       gSigmaCorrConservativePC->SetPoint(ibin,pt,0.);
 235       gSigmaCorrConservativePC->SetPointError(ibin,valFcErrx,valFcErrx,avErryfdlPC,avErryfdhPC);
 236       gFcCorrConservative->SetPoint(ibin,pt,corrfd);
 237       gFcCorrConservative->SetPointError(ibin,valFcErrx,valFcErrx,corrfdl,corrfdh);
 238     }
 239
 240   }
 241
 242
 243   gROOT->SetStyle("Plain");
 244   gStyle->SetOptTitle(0);
 245
 246   //
 247   // Plot the results
 248   TH2F *histo2Draw = new TH2F("histo2Draw","histo2 (for drawing)",100,0,20.,100,1e3,5e7);
 249   histo2Draw->SetStats(0);
 250   histo2Draw->GetXaxis()->SetTitle("p_{T}  [GeV]");
 251   histo2Draw->GetXaxis()->SetTitleSize(0.05);
 252   histo2Draw->GetXaxis()->SetTitleOffset(0.95);
 253   histo2Draw->GetYaxis()->SetTitle("#frac{1}{BR} #times #frac{d#sigma}{dp_{T}} |_{|y|<0.5}");
 254   histo2Draw->GetYaxis()->SetTitleSize(0.05);
 255   //
 256   TCanvas *combinefdunc = new TCanvas("combinefdunc","show the FD results combination");
 257   //
 258   histo2Draw->Draw();
 259   //
 260   histoSigmaCorrFc->SetMarkerStyle(20);
 261   histoSigmaCorrFc->SetMarkerColor(kGreen+2);
 262   histoSigmaCorrFc->SetLineColor(kGreen+2);
 263   histoSigmaCorrFc->Draw("esame");
 264   gSigmaCorrConservativeFc->SetMarkerStyle(20);
 265   gSigmaCorrConservativeFc->SetMarkerColor(kGreen+2);
 266   gSigmaCorrConservativeFc->SetLineColor(kGreen+2);
 267   gSigmaCorrConservativeFc->SetFillStyle(3004);//2);
 268   gSigmaCorrConservativeFc->SetFillColor(kGreen);
 269   gSigmaCorrConservativeFc->Draw("2[]same");
 270   //
 271   histoSigmaCorrNb->SetMarkerStyle(25);
 272   histoSigmaCorrNb->SetMarkerColor(kViolet+5);
 273   histoSigmaCorrNb->SetLineColor(kViolet+5);
 274   histoSigmaCorrNb->Draw("esame");
 275   gSigmaCorrConservativeNb->SetMarkerStyle(25);
 276   gSigmaCorrConservativeNb->SetMarkerColor(kOrange+7);//kViolet+5);
 277   gSigmaCorrConservativeNb->SetLineColor(kOrange+7);//kOrange+7);//kViolet+5);
 278   gSigmaCorrConservativeNb->SetFillStyle(3018);//02);
 279   gSigmaCorrConservativeNb->SetFillColor(kMagenta);
 280   gSigmaCorrConservativeNb->Draw("2[]same");
 281   //
 282   gSigmaCorrConservative->SetLineColor(kRed);
 283   gSigmaCorrConservative->SetLineWidth(2);
 284   gSigmaCorrConservative->SetFillColor(kRed);
 285   gSigmaCorrConservative->SetFillStyle(0);
 286   gSigmaCorrConservative->Draw("2");
 287   histoSigmaCorr->SetMarkerColor(kRed);
 288   histoSigmaCorr->Draw("esame");
 289   //
 290   //
 291   TLegend* leg=combinefdunc->BuildLegend();
 292   leg->SetFillStyle(0);
 293   combinefdunc->SetLogy();
 294   combinefdunc->Update();
 295
 296   TCanvas *combinefcunc = new TCanvas("combinefcunc","show the fc FD results combination");
 297   //
 298   TH2F *histo3Draw = new TH2F("histo3Draw","histo3 (for drawing)",100,0,20.,10,0.,1.);
 299   histo3Draw->SetStats(0);
 300   histo3Draw->GetXaxis()->SetTitle("p_{T}  [GeV]");
 301   histo3Draw->GetXaxis()->SetTitleSize(0.05);
 302   histo3Draw->GetXaxis()->SetTitleOffset(0.95);
 303   histo3Draw->GetYaxis()->SetTitle("Prompt fraction of the raw yields");
 304   histo3Draw->GetYaxis()->SetTitleSize(0.05);
 305   histo3Draw->Draw();
 306   //
 307   gFcConservativeFc->SetMarkerStyle(20);
 308   gFcConservativeFc->SetMarkerColor(kGreen+2);
 309   gFcConservativeFc->SetLineColor(kGreen+2);
 310   gFcConservativeFc->SetFillStyle(3004);
 311   gFcConservativeFc->SetFillColor(kGreen);
 312   gFcConservativeFc->Draw("2P");
 313   //
 314   gFcConservativeNb ->SetMarkerStyle(25);
 315   gFcConservativeNb ->SetMarkerSize(1.3);
 316   gFcConservativeNb->SetMarkerColor(kOrange+7);//kViolet+5);
 317   gFcConservativeNb->SetLineColor(kOrange+7);//kViolet+5);
 318   gFcConservativeNb->SetFillStyle(3018);
 319   gFcConservativeNb->SetFillColor(kMagenta);
 320   gFcConservativeNb->Draw("2P");
 321   //
 322   gFcCorrConservative->SetMarkerStyle(21);
 323   gFcCorrConservative->SetLineColor(kRed);
 324   gFcCorrConservative->SetLineWidth(2);
 325   gFcCorrConservative->SetFillColor(kRed);
 326   gFcCorrConservative->SetFillStyle(0);
 327   gFcCorrConservative->Draw("2P");
 328   //
 329   leg=combinefcunc->BuildLegend();
 330   leg->SetFillStyle(0);
 331   //
 332   combinefcunc->Update();
 333
 334   //
 335   // Plot the results
 336   TCanvas *finalresults = new TCanvas("finalresults","show all combined results");
 337   //
 338   if ( decay==1 ) {
 339     thD0KpifromBprediction->SetLineColor(kGreen+2);
 340     thD0KpifromBprediction->SetLineWidth(3);
 341     thD0KpifromBprediction->SetFillColor(kGreen-6);
 342     thD0KpifromBprediction->Draw("3CA");
 343     thD0KpifromBprediction->Draw("CX");
 344   }
 345   else if ( decay==2 ) {
 346     thDpluskpipiprediction->SetLineColor(kGreen+2);
 347     thDpluskpipiprediction->SetLineWidth(3);
 348     thDpluskpipiprediction->SetFillColor(kGreen-6);
 349     thDpluskpipiprediction->Draw("3CA");
 350     thDpluskpipiprediction->Draw("CX");
 351   }
 352   else if ( decay==3 ) {
 353     thDstarD0piprediction->SetLineColor(kGreen+2);
 354     thDstarD0piprediction->SetLineWidth(3);
 355     thDstarD0piprediction->SetFillColor(kGreen-6);
 356     thDstarD0piprediction->Draw("3CA");
 357     thDstarD0piprediction->Draw("CX");
 358   }
 359   else if ( decay==4 ) {
 360     thDsKKpiprediction->SetLineColor(kGreen+2);
 361     thDsKKpiprediction->SetLineWidth(3);
 362     thDsKKpiprediction->SetFillColor(kGreen-6);
 363     thDsKKpiprediction->Draw("3CA");
 364     thDsKKpiprediction->Draw("CX");
 365   }
 366   //
 367   gSigmaCorr->SetLineColor(kRed);
 368   gSigmaCorr->SetLineWidth(1);
 369   gSigmaCorr->SetFillColor(kRed);
 370   gSigmaCorr->SetFillStyle(0);
 371   gSigmaCorr->Draw("2");
 372   histoSigmaCorr->SetMarkerStyle(21);
 373   histoSigmaCorr->SetMarkerColor(kRed);
 374   histoSigmaCorr->Draw("esame");
 375   //
 376   leg = new TLegend(0.7,0.75,0.87,0.5);
 377   leg->SetBorderSize(0);
 378   leg->SetLineColor(0);
 379   leg->SetFillColor(0);
 380   leg->SetTextFont(42);
 381   if ( decay==1 ) leg->AddEntry(thD0KpifromBprediction,"FONLL ","fl");
 382   else if ( decay==2 ) leg->AddEntry(thDpluskpipiprediction,"FONLL ","fl");
 383   else if ( decay==3 ) leg->AddEntry(thDstarD0piprediction,"FONLL ","fl");
 384   else if ( decay==4 ) leg->AddEntry(thDsKKpiprediction,"FONLL ","fl");
 385   leg->AddEntry(histoSigmaCorr,"data stat. unc.","pl");
 386   leg->AddEntry(gSigmaCorr,"data syst. unc.","f");
 387   leg->Draw();
 388   //
 389   finalresults->SetLogy();
 390   finalresults->Update();
 391
 392
 393   //
 394   // Draw all the systematics independently
 395   systematics.DrawErrors(gSigmaCorrConservativePC);
 396
 397
 398   // Write the output to a file
 399   TFile * out = new TFile(outfilename,"recreate");
 400   histoSigmaCorr->Write();
 401   gSigmaCorr->Write();
 402   gSigmaCorrConservative->Write();
 403   gSigmaCorrConservativePC->Write();
 404   gFcCorrConservative->Write();
 405   out->Write();
 406
 407 }