PWG4/macros/jetspectrum/runUnfoldingCF.C

   1 void Load();
   2 TList *GetResults(Char_t *testfile,char* listName);
   3 TObject* GetMeasuredSpectrum(TList *fContainer);
   4 TObject* GetGeneratedSpectrum(TList *fContainer);
   5 TObject* GetEfficiency(TList *fContainer,Int_t iCase = 44);
   6 THnSparse* GetResponseMatrix(TList *fContainer);
   7 THnSparse* CreateResponseMatrix();
   8 Float_t GetNTrials(TList *fContainer);
   9 THnSparse* CreateGuessed(const THnSparse* h);
  10 void PrintErrors(THnSparse *h);
  11
  12 Int_t fDebug = 10;
  13
  14 int iRebin = 8; //
  15
  16 // *************************************************************************************************************
  17 // Macro for unfolding jet spectra:
  18 // * First Unfolding is always the MC itself
  19 // * if bMCtest we take a second MC which is unfolded with the first MC otherwise real data is unfolded with first MC
  20 // * bMC2 controls wether we correct to charged particle jets or not
  21 // ############### CHANGELOG ################
  22 // 23.11. added method to create an empty response matrix
  23
  24 void runUnfoldingCF( int gIterations = 10, bool bMC2 = true,bool bNoPrior = true,bool bMCtest = true, int iJetF = 0) {
  25   Load();
  26
  27   const bool bPreSmooth = false;
  28
  29   // can be different depeding on the prior
  30   Int_t fIterations = gIterations;
  31   Int_t fIterations2 = gIterations;
  32   Int_t fIterations3 = gIterations;
  33
  34   // the first one is usually MC only
  35   TString cJetF;
  36   TString listName;
  37   TString testfile;
  38   if(bMCtest) testfile = "allpt_lhc10e14_100903_Chunk1.root"; // the higher stats chunck
  39   else testfile = "allpt_lhc10e14_100903.root";
  40   TList *results = 0;
  41   if(iJetF == 0){
  42     cJetF = "anti-k_{T}";
  43     listName = "spec2_jetsAOD_FASTJET04_jetsAODMC%s_FASTJET04_0000000000";
  44   }
  45   else if(iJetF == 1){
  46     cJetF = "k_{T}";
  47     listName = "spec2_jetsAOD_FASTKT04_jetsAODMC%s_FASTKT04_0000000000";
  48   }
  49   else if(iJetF == 2){
  50     cJetF = "UA1";
  51     listName = "spec2_jets_jetsAODMC%s_UA104_0000000000";
  52   }
  53   TList *results =  GetResults(testfile.Data(),Form(listName.Data(),(bMC2?"2":"")));
  54
  55   // here comes the real data
  56   TString cJetF2;
  57   TString listName2;
  58   TString testfile2;
  59
  60   if(bMCtest)testfile2 = "allpt_lhc10e14_100903_Chunk0.root";  // the smalle stats
  61   else testfile2 = "/Users/kleinb/alice/jets/train/100910/PWG4_JetTasksOutput_Merge_bcd.root";
  62   if(iJetF == 0){
  63     cJetF = "anti-k_{T}";
  64     if(bMCtest)listName2 = "spec2_jetsAOD_FASTJET04_jetsAODMC%s_FASTJET04_0000000000";
  65     else listName2 = "spec2_jetsAOD_FASTJET04__0000000000";
  66
  67   }
  68   else if(iJetF == 1){
  69     cJetF = "k_{T}";
  70     if(bMCtest)listName2 = "spec2_jetsAOD_FASTKT04_jetsAODMC%s_FASTKT04_0000000000";
  71     else listName2 = "spec2_jetsAOD_FASTKT04__0000000000";
  72
  73   }
  74   else if(iJetF == 2){
  75     cJetF = "k_{T}";
  76     if(bMCtest)listName2 = "spec2_jets_jetsAODMC%s_UA104_0000000000";
  77     else listName2 = "spec2_jets__0000000000";
  78
  79   }
  80   TList *results2 = GetResults(testfile2.Data(),Form(listName2.Data(),(bMC2?"2":"")));bool bScale = false;
  81
  82
  83
  84   if(!results){
  85     Error("No output objects: Calculation will terminate here");
  86     return;
  87   }
  88
  89   bool b1D = true;
  90   const Float_t kMaxPt = 140; //  only for drawin...
  91   const Float_t kMaxPt2 = 140;
  92
  93
  94   char cString[] = Form("unfolded/101010_unfolding%s_MC%s%s%%s_%sIter%04d_Rebin%02d_jet%d",(bNoPrior?"_NoPrior":""),(bMC2?"2":""),(bPreSmooth?"_PreSmooth":""),(bMCtest?"MCtest_":""),gIterations,iRebin,iJetF);
  95   char cPrintMask[] = Form("%s.png",cString);
  96
  97   TF3 *fSmooth = new TF3("fSmooth","[0]*pow(x,[1])+[2]*y+[3]*z",5.,150,-10,10,-10,10);
  98   fSmooth = 0;
  99   TF3 *fSmooth2 = 0;
 100   TF3 *fSmooth3 = 0;
 101
 102   //  fSmooth = 0;
 103   AliLog::SetGlobalLogLevel(AliLog::kDebug);
 104   // get the essential
 105   if(fSmooth){
 106     fSmooth->SetParameters(1,-8,0,0);
 107     fSmooth2 = (TF3*)fSmooth->Clone("fSmooth2");
 108     fSmooth3 = (TF3*)fSmooth->Clone("fSmooth3");
 109   }
 110
 111
 112
 113
 114   THnSparseF *efficiency  = (THnSparseF*)  GetEfficiency(results,41);
 115   THnSparseF *efficiencyRecMatch  = (THnSparseF*)  GetEfficiency(results,131);
 116   THnSparseF *response    = (THnSparseF*)  GetResponseMatrix(results);
 117
 118
 119
 120   THnSparseF *measuredIn2    = (THnSparseF*) GetMeasuredSpectrum(results2);
 121   measuredIn2->Multiply(efficiencyRecMatch);
 122   THnSparseF *measuredIn    = (THnSparseF*) GetMeasuredSpectrum(results);
 123   measuredIn->Multiply(efficiencyRecMatch);
 124   THnSparseF *generatedOut   = (THnSparseF*) GetGeneratedSpectrum(results);
 125   THnSparseF *generatedOut2   = (THnSparseF*) GetGeneratedSpectrum(results2);
 126
 127   Float_t fTrials =  GetNTrials(results);
 128   Float_t fTrials2 = GetNTrials(results2);
 129   if(bScale){
 130     Printf("Trials %.3E %.3E",fTrials,fTrials2);
 131     measuredIn2->Scale(1./fTrials2);
 132     measuredIn->Scale(1./fTrials);
 133   }
 134
 135   // set the content to zero above threshold
 136   /*
 137   const float lastPt = 250;
 138   Int_t ibxyz[3];
 139   for(int ibx = 1;ibx<=measuredIn2->GetAxis(0)->GetNbins();ibx++){
 140     float pt = measuredIn2->GetAxis(0)->GetBinCenter(ibx);
 141
 142     if(pt<lastPt)continue;
 143     ibxyz[0] = ibx;
 144     for(int iby = 1;iby<=measuredIn2->GetAxis(1)->GetNbins();iby++){
 145       ibxyz[1] = iby;
 146       for(int ibz = 1;ibz<=measuredIn2->GetAxis(2)->GetNbins();ibz++){
 147         ibxyz[2] = ibz;
 148         measuredIn2->SetBinContent(ibxyz,0);
 149         measuredIn2->SetBinError(ibxyz,0);
 150       }
 151     }
 152   }
 153   */
 154
 155
 156   // for testing
 157   const int nDim = 3;
 158   const int dimrec[nDim] = {0,1,2};
 159   const int dimgen[nDim] = {3,4,5};
 160
 161
 162   /*
 163   THnSparseF *generated = response->Projection(nDim,dimgen);
 164   THnSparseF *measured = response->Projection(nDim,dimrec);
 165   */
 166
 167   // create a guessed "a priori" distribution using binning of MC
 168   THnSparse* guessed2 = CreateGuessed(generatedOut) ; // can at best take the measured?
 169   THnSparse* guessed = CreateGuessed(generatedOut) ; // can at best take the measured?
 170   Printf("%s:%d %d",(char*)__FILE__,__LINE__,guessed2->GetNbins());
 171   Printf("%s:%d %d",(char*)__FILE__,__LINE__,guessed->GetNbins());
 172
 173   //---- Dbug show the errrors
 174   //  PrintErrors(guessed);
 175   //  PrintErrors(response);
 176   //  PrintErrors(efficiency);
 177   //   PrintErrors(measuredIn);
 178   //  PrintErrors(generatedOut);
 179
 180   Bool_t bCorrelatedErrors = true;
 181   AliCFUnfolding unfolding("unfolding","title",3,response,efficiency,measuredIn,(bNoPrior?0:guessed));
 182   unfolding.SetMaxNumberOfIterations(fIterations); // regulate flutuations...
 183   unfolding.SetMaxChi2PerDOF(0);
 184   if(fSmooth)unfolding.UseSmoothing(fSmooth);
 185   if(bCorrelatedErrors){
 186     unfolding.SetUseCorrelatedErrors(kTRUE);
 187     unfolding.SetMaxConvergencePerDOF(0.01*0.01);
 188   }
 189
 190   unfolding.Unfold();
 191
 192   THnSparse* result = unfolding.GetUnfolded();
 193   THnSparse* estMeasured = unfolding.GetEstMeasured();
 194
 195
 196   //----
 197   AliCFUnfolding unfolding2("unfolding2","title",3,response,efficiency,measuredIn2,(bNoPrior?0:guessed));
 198   //  AliCFUnfolding unfolding2("unfolding2","title",3,response,efficiency,measuredIn2,0);
 199   unfolding2.SetMaxNumberOfIterations(fIterations2);
 200   unfolding2.SetMaxChi2PerDOF(0);
 201   // carefull with 0x0 pointer neighbouring bin smoothing is switched on...
 202   if(fSmooth2)unfolding2.UseSmoothing(fSmooth2);
 203   if(bCorrelatedErrors){
 204     unfolding2.SetUseCorrelatedErrors(kTRUE);
 205     unfolding2.SetMaxConvergencePerDOF(0.01*0.01);
 206   }
 207   unfolding2.Unfold();
 208
 209
 210
 211   THnSparse* estMeasuredTmp = unfolding2.GetEstMeasured();
 212   THnSparse* result2 = 0;
 213   THnSparse* estMeasured2 = 0;
 214   if(bPreSmooth){
 215     // use the result of the previous as smoothed input...
 216     AliCFUnfolding unfolding3("unfolding3","title",3,response,efficiency,estMeasuredTmp,(bNoPrior?0:guessed));
 217     unfolding3.SetMaxNumberOfIterations(fIterations3);
 218     unfolding3.SetMaxChi2PerDOF(0);
 219     unfolding3.UseSmoothing(fSmooth3);
 220     unfolding3.Unfold();
 221     result2 = unfolding3.GetUnfolded();
 222     estMeasured2 = unfolding3.GetEstMeasured();
 223   }
 224   else{
 225     result2 = unfolding2.GetUnfolded();
 226     estMeasured2 = estMeasuredTmp;
 227   }
 228
 229
 230
 231   TCanvas * canvas = new TCanvas("canvas","",1200,900);
 232   canvas->Divide(2,3);
 233   TCanvas * cPrint = new TCanvas("cPrint","");
 234
 235
 236
 237   TFile *f1 = new TFile(Form("%s.root",Form(cString,"")),"RECREATE");
 238   TParameter<Long_t>* fIterationsPara =  new TParameter<Long_t> ("fIterations", 0);
 239   fIterationsPara->SetVal((Long_t)gIterations);
 240   fIterationsPara->Write();
 241
 242   TDirectory *dMC = f1->mkdir("unfoldMC");
 243   TDirectory *dReal = f1->mkdir("unfoldReal");
 244
 245   gROOT->cd();
 246   // color code black is unfolded
 247   // blue is measured
 248   // red  is generated
 249   // kGreen is guessed
 250
 251
 252
 253   if(b1D){
 254
 255     TH1* h_gen = generatedOut->Projection(0);
 256
 257     h_gen->SetMarkerColor(kRed);
 258     h_gen->SetMarkerStyle(kFullSquare);
 259     h_gen->SetName("generated");
 260     h_gen->SetTitle("generated");
 261     TH1* h_meas = measuredIn->Projection(0);
 262     h_meas->SetMarkerColor(kBlue);
 263     h_meas->SetMarkerStyle(kFullSquare);
 264     h_meas->SetName("measured");
 265     h_meas->SetTitle("measured");
 266     TH1* h_guess = guessed->Projection(0);
 267     h_guess->SetMarkerColor(kGreen);
 268     h_guess->SetMarkerStyle(kFullSquare);
 269     h_guess->SetName("guessed");
 270     h_guess->SetTitle("guesse");
 271     TH1* h_unf = result->Projection(0);
 272     h_unf->SetMarkerColor(kBlack);
 273     h_unf->SetMarkerStyle(kFullSquare);
 274     h_unf->SetName("unfolded");
 275     h_unf->SetTitle("unfolded");
 276     TH1* h_estmeas = estMeasured->Projection(0);
 277     h_estmeas->SetMarkerColor(kGray);
 278     h_estmeas->SetMarkerStyle(kFullSquare);
 279     h_estmeas->SetName("estmeas");
 280     h_estmeas->SetTitle("estmeas");
 281
 282     TLegend *leg = new TLegend(0.6,0.6,0.85,0.8);
 283     leg->SetFillColor(0);
 284     leg->SetTextFont(gStyle->GetTextFont());
 285     leg->SetBorderSize(0);
 286     leg->AddEntry(h_meas,"measured","P");
 287     leg->AddEntry(h_gen,"generated","P");
 288     leg->AddEntry(h_unf,"unfolded","P");
 289     leg->AddEntry(h_estmeas,"R * u","P");
 290
 291
 292
 293     canvas->cd(1)->SetLogy();
 294     h_gen->SetAxisRange(0,kMaxPt);
 295     h_meas->SetAxisRange(0,kMaxPt);
 296     h_meas->SetXTitle("p_{T} (GeV)");
 297     h_meas->SetYTitle("yield (arb. units)");
 298     h_meas->DrawCopy("P");
 299     h_gen->DrawCopy("Psame");
 300     h_unf->DrawCopy("Psame");
 301     h_estmeas->DrawCopy("Psame");
 302     leg->Draw("");
 303     cPrint->cd()->SetLogy();
 304     h_meas->DrawCopy("P");
 305     h_gen->DrawCopy("Psame");
 306     h_unf->DrawCopy("Psame");
 307     h_estmeas->DrawCopy("Psame");
 308     leg->Draw("");
 309     canvas->Update();
 310     cPrint->Update();
 311     cPrint->SaveAs(Form(cPrintMask,"spectrum_mc"));
 312     if(!gROOT->IsBatch()){
 313       if(getchar()=='q')return;
 314     }
 315     // the same for the real data
 316     TH1* h_gen2 = generatedOut2->Projection(0);
 317
 318     h_gen2->SetMarkerColor(kRed);
 319     h_gen2->SetMarkerStyle(kFullCircle);
 320     h_gen2->SetName("generated2");
 321     h_gen2->SetTitle("generated2");
 322     TH1* h_meas2 = measuredIn2->Projection(0);
 323     h_meas2->SetMarkerColor(kBlue);
 324     h_meas2->SetMarkerStyle(kFullCircle);
 325     h_meas2->SetName("measured2");
 326     h_meas2->SetTitle("measured2");
 327     TH1* h_guess2 = guessed2->Projection(0);
 328     h_guess2->SetMarkerColor(kGreen);
 329     h_guess2->SetMarkerStyle(kFullCircle);
 330     h_guess2->SetName("guessed2");
 331     h_guess2->SetTitle("guesse2");
 332     TH1* h_unf2 = result2->Projection(0);
 333     PrintErrors(result2);
 334     h_unf2->SetMarkerColor(kBlack);
 335     h_unf2->SetMarkerStyle(kFullCircle);
 336     h_unf2->SetName("unfolded2");
 337     h_unf2->SetTitle("unfolded2");
 338     TH1* h_estmeas2 = estMeasured2->Projection(0);
 339     h_estmeas2->SetMarkerColor(kGray);
 340     h_estmeas2->SetMarkerStyle(kFullCircle);
 341     h_estmeas2->SetName("estmeas2");
 342     h_estmeas2->SetTitle("estmeas2");
 343
 344     TLegend *leg2 = new TLegend(0.6,0.6,0.85,0.8);
 345     leg2->SetFillColor(0);
 346     leg2->SetTextFont(gStyle->GetTextFont());
 347     leg2->SetBorderSize(0);
 348     leg2->AddEntry(h_meas2,"measured","P");
 349     leg2->AddEntry(h_gen2,"generated","P");
 350     leg2->AddEntry(h_unf2,"unfolded","P");
 351     leg2->AddEntry(h_estmeas2,"R * u","P");
 352
 353     Printf("Integral Measured: %E Unfolded %E R*U %E",h_meas2->Integral(),h_unf2->Integral(),h_estmeas2->Integral());
 354
 355
 356
 357     canvas->cd(2)->SetLogy();
 358     //
 359
 360     h_meas2->SetXTitle("p_{T} (GeV)");
 361     h_meas2->SetYTitle("yield (arb. units)");
 362     //    h_gen2->DrawCopy("P");
 363     h_unf2->SetAxisRange(0,kMaxPt2);
 364     h_meas2->SetAxisRange(0,kMaxPt2);
 365     h_meas2->DrawCopy("P");
 366     h_unf2->DrawCopy("Psame");
 367     h_estmeas2->DrawCopy("Psame");
 368     leg2->Draw();
 369     cPrint->cd()->SetLogy();
 370     h_unf2->SetAxisRange(0,kMaxPt2);
 371     h_meas2->SetAxisRange(0,kMaxPt2);
 372     h_meas2->DrawCopy("P");
 373     h_unf2->DrawCopy("Psame");
 374     h_estmeas2->DrawCopy("Psame");
 375     canvas->Update();
 376     cPrint->Update();
 377     leg2->Draw();
 378     cPrint->SaveAs(Form(cPrintMask,"spectrum_real"));
 379     if(!gROOT->IsBatch()){
 380       if(getchar()=='q')return;
 381     }
 382
 383
 384
 385     // residuals
 386
 387     TH1D *hResiduals = (TH1D*) h_meas->Clone("hResiduals");
 388     hResiduals->Reset();
 389     for(int ib = 1;ib <  hResiduals->GetNbinsX();ib++){
 390       float val1 =  h_meas->GetBinContent(ib);
 391       float err1 =  h_meas->GetBinError(ib);
 392       float val2 =  h_estmeas->GetBinContent(ib);
 393       if(err1>0){
 394         Float_t res = (val1-val2)/err1;
 395         Float_t res_err = (val1-val2)/err1*0.01; // error bars of 1%
 396         hResiduals->SetBinContent(ib,res);
 397         hResiduals->SetBinError(ib,0.01);
 398       }
 399
 400     }
 401     hResiduals->SetXTitle("p_{T} (GeV)");
 402     hResiduals->SetYTitle("#frac{m - R * u}{e_{m}}");
 403     hResiduals->SetMaximum(4.5);
 404     hResiduals->SetMinimum(-4.5);
 405
 406
 407     canvas->cd(3);
 408     //    hRatioGenUnf->DrawCopy("p");
 409     hResiduals->SetAxisRange(0,kMaxPt);
 410     hResiduals->DrawCopy("P");
 411     cPrint->cd();
 412     gPad->SetLogy(0);
 413     hResiduals->DrawCopy("P");
 414     canvas->Update();
 415     cPrint->Update();
 416     cPrint->SaveAs(Form(cPrintMask,"residuals_mc"));
 417     if(!gROOT->IsBatch()){
 418       if(getchar()=='q')return;
 419     }
 420
 421     TH1D *hResiduals2 =  (TH1D*)h_meas2->Clone("hResiduals2");
 422     hResiduals2->Reset();
 423     for(int ib = 1;ib <  hResiduals2->GetNbinsX();ib++){
 424       float val1 =  h_meas2->GetBinContent(ib);
 425       float err1 =  h_meas2->GetBinError(ib);
 426       float val2 =  h_estmeas2->GetBinContent(ib);
 427       if(err1>0){
 428         Float_t res = (val1-val2)/err1;
 429         Float_t res_err = (val1-val2)/err1*0.01; // error bars of 1%
 430         hResiduals2->SetBinContent(ib,res);
 431         hResiduals2->SetBinError(ib,0.01);
 432       }
 433
 434     }
 435     hResiduals2->SetMaximum(4.5);
 436     hResiduals2->SetMinimum(-4.5);
 437     hResiduals2->SetYTitle("#frac{m - R * u}{e_{m}}");
 438     hResiduals2->SetAxisRange(0,kMaxPt2);
 439     cPrint->cd();
 440     gPad->SetLogy(0);
 441     hResiduals2->DrawCopy("P");
 442     canvas->Update();
 443     cPrint->Update();
 444     cPrint->SaveAs(Form(cPrintMask,"residuals_real"));
 445     if(!gROOT->IsBatch()){
 446       if(getchar()=='q')return;
 447     }
 448
 449     canvas->Update();
 450     cPrint->Update();
 451
 452     if(!gROOT->IsBatch()){
 453       if(getchar()=='q')return;
 454     }
 455     hResiduals2->DrawCopy("P");
 456     canvas->Update();
 457     cPrint->Update();
 458     if(!gROOT->IsBatch()){
 459       if(getchar()=='q')return;
 460     }
 461
 462     TH1D *hRatioGenUnf2 = (TH1D*) h_gen2->Clone("hRatioGenUnf2");
 463     hRatioGenUnf2->Divide(h_unf2);
 464     canvas->cd(4);
 465     //    hRatioGenUnf2->DrawCopy("p");
 466
 467     hResiduals2->DrawCopy("P");
 468
 469     TH1D *hRatioUnfMeas = (TH1D*) h_unf->Clone("hRatioUnfMeas");
 470     hRatioUnfMeas->SetXTitle("p_{T} (GeV)");
 471     hRatioUnfMeas->SetYTitle("unfolded/meas");
 472     hRatioUnfMeas->Divide(h_meas);
 473     hRatioUnfMeas->SetMaximum(15);
 474     if(bMC2)hRatioUnfMeas->SetMaximum(5);
 475     hRatioUnfMeas->SetMinimum(0);
 476     TH1D *hRatioGenMeas =  (TH1D*)h_gen->Clone("hRatioGenMeas");
 477     hRatioGenMeas->Divide(h_meas);
 478
 479     canvas->cd(5);
 480
 481     hRatioUnfMeas->SetAxisRange(0.,kMaxPt);
 482     hRatioUnfMeas->DrawCopy("p");
 483     hRatioGenMeas->DrawCopy("psame");
 484
 485
 486
 487
 488     TH1D *hRatioUnfMeas2 =  (TH1D*) h_unf2->Clone("hRatioUnfMeas2");
 489     hRatioUnfMeas2->Divide(h_meas2);
 490     canvas->cd(6);
 491     hRatioUnfMeas2->SetXTitle("p_{T} (GeV)");
 492     hRatioUnfMeas2->SetYTitle("unfolded/meas");
 493     hRatioUnfMeas2->SetAxisRange(0.,kMaxPt2);
 494     hRatioUnfMeas2->SetMaximum(15);
 495     if(bMC2)hRatioUnfMeas2->SetMaximum(5);
 496     hRatioUnfMeas2->SetMinimum(0);
 497     hRatioUnfMeas2->DrawCopy("p");
 498     hRatioUnfMeas->DrawCopy("psame");
 499
 500     // plot the effieciencies
 501     cPrint->cd();
 502     TH1* hEffGen = efficiency->Projection(0);
 503     hEffGen->SetName("hEffGen");
 504     TH1* hEffRec = efficiencyRecMatch->Projection(0);
 505     hEffRec->SetName("hEffRec");
 506
 507     hEffGen->SetAxisRange(0,kMaxPt2);
 508     hEffRec->SetAxisRange(0,kMaxPt2);
 509
 510     hEffRec->SetXTitle("p_{T,gen/rec} (GeV/c)");
 511     hEffRec->SetYTitle("Efficiency");
 512     hEffGen->SetXTitle("p_{T,gen/rec} (GeV/c)");
 513     hEffGen->SetYTitle("Efficiency");
 514
 515     hEffGen->SetMarkerColor(kBlue);
 516     hEffGen->SetMarkerStyle(kFullCircle);
 517     hEffRec->SetMarkerColor(kBlue);
 518     hEffRec->SetMarkerStyle(kOpenCircle);
 519     hEffGen->SetMaximum(1.2);
 520     hEffGen->SetMinimum(0.3);
 521     hEffGen->DrawCopy("p");
 522     hEffRec->DrawCopy("psame");
 523     cPrint->Update();
 524     cPrint->SaveAs(Form(cPrintMask,"effs"));
 525
 526     h_unf->SetDirectory(dMC);
 527     h_meas->SetDirectory(dMC);
 528     hResiduals->SetDirectory(dMC);
 529     h_estmeas->SetDirectory(dMC);
 530     h_gen->SetDirectory(dMC);
 531
 532
 533     // put the used effs and response to the unfolded of the MC
 534
 535     hEffGen->SetDirectory(dMC);
 536     hEffRec->SetDirectory(dMC);
 537     //    response->SetDirectory(dMC);
 538     dMC->cd();
 539     guessed->Write();
 540     response->Write();
 541     measuredIn->Write();
 542     gROOT->cd();
 543
 544
 545     h_unf2->SetDirectory(dReal);
 546     h_meas2->SetDirectory(dReal);
 547     hResiduals2->SetDirectory(dReal);
 548     h_estmeas2->SetDirectory(dReal);
 549     h_gen2->SetDirectory(dReal); // will be empty if we use real data, but can check with different data set as well...
 550
 551     dReal->cd();
 552     guessed2->Write();
 553     measuredIn2->Write();
 554     gROOT->cd();
 555     // store the parameters of the unfloding with tparamter...
 556
 557
 558   }
 559   canvas->cd();
 560   canvas->Modified();
 561   canvas->SaveAs(Form(cPrintMask,"all"));
 562
 563   f1->Write();
 564   f1->Close();
 565
 566
 567   return;
 568   if(fDebug)Printf("Line:%d %s Entries %f",__LINE__,h_gen->GetName(),h_gen->GetEntries());
 569   //printf("c1\n");
 570   if(b1D)h_gen->DrawCopy("P");
 571   else h_gen->DrawCopy("lego2");
 572
 573
 574
 575
 576   if(b1D) canvas->cd(2)->SetLogy();
 577   else canvas->cd(2)->SetLogz();
 578   TH1* h_meas = 0;
 579   if(b1D) h_meas = measuredIn->Projection(0);
 580   else h_guessed = measuredIn->Projection(0,1);
 581   h_meas->SetName("measured");
 582   h_meas->SetTitle("measured");
 583   if(fDebug)Printf("Line:%d %s Entries %f",__LINE__,h_meas->GetName(),h_meas->GetEntries());
 584
 585   //printf("c2\n");
 586   if(b1D) h_meas->Draw("hist");
 587   else h_meas->Draw("lego2");
 588
 589
 590   if(b1D) canvas->cd(3)->SetLogy();
 591   else canvas->cd(3)->SetLogz();
 592
 593   TH1* h_guessed = 0;
 594   if(b1D) h_guessed = guessed2->Projection(0);
 595   else h_guessed = guessed2->Projection(0,1);
 596
 597   h_guessed->SetTitle("a priori");
 598   h_guessed->SetName("apriori");
 599   if(fDebug)Printf("Line:%d %s Entries %f",__LINE__,h_guessed->GetName(),h_guessed->GetEntries());
 600   //printf("c3\n");
 601   if(b1D) h_guessed->Draw("hist");
 602   else h_guessed->Draw("lego2");
 603
 604   canvas->cd(4);
 605   TH1* h_eff = 0;
 606   if(b1D) h_eff = efficiency->Projection(0);
 607   else h_eff = efficiency->Projection(0,1);
 608
 609   h_eff->SetTitle("efficiency");
 610   h_eff->SetName("efficiency");
 611   if(fDebug)Printf("Line:%d %s Entries %f",__LINE__,h_eff->GetName(),h_eff->GetEntries());
 612   //printf("c4\n");
 613   if(b1D) h_eff->Draw("hist");
 614   else h_eff->Draw("lego2");
 615
 616
 617   if(b1D) canvas->cd(5)->SetLogy();
 618   else canvas->cd(5)->SetLogz();
 619   TH1* h_unf = 0;
 620   if(b1D) h_unf = result->Projection(0);
 621   else h_unf = result->Projection(0,1);
 622
 623   h_unf->SetName("unfolded");
 624   h_unf->SetTitle("unfolded");
 625   //printf("c5\n");
 626   if(b1D) h_unf->Draw("hist");
 627   else h_unf->Draw("lego2");
 628
 629   canvas->cd(6);
 630   TH1* ratio = 0;
 631   if(b1D){
 632     ratio = (TH1D*)h_unf->Clone("ratio");
 633     ratio->SetTitle("unfolded / generated");
 634     ratio->Divide(h_unf,h_gen,1,1);
 635   }
 636   else{
 637     ratio = (TH2D*)h_unf->Clone("ratio");
 638     ratio->SetTitle("unfolded / generated");
 639     ratio->Divide(h_unf,h_gen,1,1);
 640   }
 641 //   ratio->Draw("cont4z");
 642 //   ratio->Draw("surf2");
 643   //printf("c6\n");
 644   if(b1D)ratio->Draw("hist");
 645   else ratio->Draw("colz");
 646
 647   return;
 648
 649   canvas->cd(7);
 650   TH2* orig = unfolding.GetOriginalPrior()->Projection(1,0);
 651   orig->SetName("originalprior");
 652   orig->SetTitle("original prior");
 653   //printf("c7\n");
 654   orig->Draw("lego2");
 655
 656   canvas->cd(8);
 657   THnSparseF* corrected = (THnSparseF*)measured->Clone("corrected");
 658   //corrected->ApplyEffCorrection(*efficiency);
 659   TH2D* corr = corrected->Projection(0,1);
 660   corr->SetTitle("simple correction");
 661   corr->SetName("simplecorrection");
 662   //printf("c8\n");
 663   corr->Draw("lego2");
 664
 665   canvas->cd(9);
 666   TH2D* ratio2 = (TH2D*) corr->Clone();
 667   ratio2->Divide(corr,h_gen,1,1);
 668   ratio2->SetTitle("simple correction / generated");
 669   //printf("c9\n");
 670   ratio2->Draw("lego2");
 671
 672   return;
 673 }
 674
 675 // ====================================================================
 676
 677
 678 void Load(){
 679   gSystem->Load("libTree");
 680   gSystem->Load("libPhysics");
 681   gSystem->Load("libHist");
 682   gSystem->Load("libVMC");
 683   gSystem->Load("libSTEERBase");
 684   gSystem->Load("libAOD");
 685   gSystem->Load("libESD");
 686   gSystem->Load("libANALYSIS");
 687   gSystem->Load("libANALYSISalice");
 688   gSystem->Load("libJETAN");
 689   gSystem->Load("libCORRFW");
 690   gSystem->Load("libPWG4JetTasks");
 691   gStyle->SetPalette(1);
 692   gStyle->SetOptStat(0);
 693 }
 694
 695 TList *GetResults( Char_t *testfile, Char_t *cName){
 696   //
 697   // read output
 698   //
 699   TFile *f = TFile::Open(testfile);
 700   if(!f || f->IsZombie()){
 701     Error("File not readable");
 702     return 0x0;
 703   }
 704   TDirectory *dir = dynamic_cast<TDirectory*>(f->Get(Form("PWG4_%s",cName)));
 705   if(!dir){
 706     Printf("%s:%d: Output directory %s not found",(char*)__FILE__,__LINE__,Form("PWG4_%s",cName));
 707     f->ls();
 708     f->Close(); delete f;
 709     return 0x0;
 710   }
 711   TList *l = dynamic_cast<TList *>(dir->Get(Form("pwg4%s",cName)));
 712   if(!l){
 713     Printf("%s:%d: Output list %s not found",(char*)__FILE__,__LINE__,Form("pwg4%s",cName));
 714     dir->ls();
 715     f->Close(); delete f;
 716     return 0x0;
 717   }
 718   //  TList *returnlist = dynamic_cast<TList *>(l->Clone());
 719   // f->Close(); delete f;
 720   return l;
 721 }
 722
 723
 724 TObject* GetMeasuredSpectrum(TList *fContainer) {
 725   THnSparseF* htmp = (THnSparseF*)fContainer->FindObject(Form("fhnJetContainer%d",AliAnalysisTaskJetSpectrum2::kMaxStep+AliAnalysisTaskJetSpectrum2::kStep1));
 726   THnSparseF* hm = (THnSparseF*)htmp->Clone("measured");
 727   hm->SetName("measured");
 728   hm->SetTitle("measured");
 729   if(iRebin==1)return hm;
 730   Int_t fRebin[3];
 731   fRebin[0] = iRebin;
 732   fRebin[1] = 1;
 733   fRebin[2] = 1;
 734   THnSparse *hm2 = hm->Rebin(fRebin); // this creates a new histo...
 735   //  hm2->Scale(iRebin);
 736   if(fDebug)Printf("Line: %d %s Entries %f",__LINE__,hm->GetName(),hm->GetEntries());
 737   return hm2;
 738 }
 739
 740 Float_t GetNTrials(TList *fContainer){
 741   TH1F* htmp = (TH1F*)fContainer->FindObject("fh1Trials");
 742   return htmp->Integral();
 743 }
 744
 745 TObject* GetGeneratedSpectrum(TList *fContainer) {
 746   // the generated jet spectrum within eta < 0.5
 747   THnSparseF* htmp = (THnSparseF*)fContainer->FindObject(Form("fhnJetContainer%d",AliAnalysisTaskJetSpectrum2::kStep1));
 748   THnSparseF* hg = (THnSparseF*)htmp->Clone("generated");
 749   hg->SetName("generated");
 750   hg->SetTitle("generated");
 751
 752   if(iRebin==1)return hg;
 753   Int_t fRebin[3];
 754   fRebin[0] = iRebin;
 755   fRebin[1] = 1;
 756   fRebin[2] = 1;
 757   THnSparse* hg2 = hg->Rebin(fRebin);
 758   // hg2->Scale(iRebin);
 759   if(fDebug)Printf("Line: %d %s Entries %f",__LINE__,hg->GetName(),hg->GetEntries());
 760   return hg2;
 761 }
 762
 763 TObject* GetEfficiency(TList *fContainer,Int_t iCase) {
 764
 765   static int count = 0;
 766   THnSparseF* hn1 = 0;
 767   THnSparseF* hn2 = 0;
 768   THnSparseF* he = 0;
 769   if(iCase==44){
 770     // Unitiy efficiency...
 771     hn1 = (THnSparseF*)fContainer->FindObject(Form("fhnJetContainer%d",AliAnalysisTaskJetSpectrum2::kStep4));
 772     hn2 = (THnSparseF*)fContainer->FindObject(Form("fhnJetContainer%d",AliAnalysisTaskJetSpectrum2::kStep4));
 773    }
 774   else if(iCase==14){
 775     // eficiency Generated with gen < 0.5 --> gen with Rec partner in 0.5
 776     hn1 = (THnSparseF*)fContainer->FindObject(Form("fhnJetContainer%d",AliAnalysisTaskJetSpectrum2::kStep1));
 777      hn2 = (THnSparseF*)fContainer->FindObject(Form("fhnJetContainer%d",AliAnalysisTaskJetSpectrum2::kStep4));
 778    }
 779   else if(iCase==41){
 780     // eficiency Generated with gen < 0.5 --> gen with Rec partner in 0.5
 781     hn1 = (THnSparseF*)fContainer->FindObject(Form("fhnJetContainer%d",AliAnalysisTaskJetSpectrum2::kStep4));
 782     hn2 = (THnSparseF*)fContainer->FindObject(Form("fhnJetContainer%d",AliAnalysisTaskJetSpectrum2::kStep1));
 783    }
 784   else if(iCase==131){
 785     // eficiency reconstruted with < 0.5 --> with partner
 786     hn1 = (THnSparseF*)fContainer->FindObject(Form("fhnJetContainer%d",AliAnalysisTaskJetSpectrum2::kMaxStep+AliAnalysisTaskJetSpectrum2::kStep3));
 787    hn2 = (THnSparseF*)fContainer->FindObject(Form("fhnJetContainer%d",AliAnalysisTaskJetSpectrum2::kMaxStep+AliAnalysisTaskJetSpectrum2::kStep1));
 788    }
 789
 790   if(iRebin==1){
 791     he = (THnSparseF*)hn1->Clone("efficiency");
 792     he->Divide(hn1,hn2,1.,1.,"B");
 793     he->SetName(Form("efficiency_%d_%d",iCase,count));
 794     he->SetTitle(Form("efficiency_%d_%d",iCase,count));
 795     count++;
 796     return he;
 797   }
 798   Int_t fRebin[3];
 799   fRebin[0] = iRebin;
 800   fRebin[1] = 1;
 801   fRebin[2] = 1;
 802
 803   THnSparseF* hn1_2 = hn1->Rebin(fRebin);
 804   THnSparseF* hn2_2 = hn2->Rebin(fRebin);
 805   he = (THnSparseF*)hn1_2->Clone("efficiency");
 806   he->Divide(hn1_2,hn2_2,1.,1.,"B");
 807   he->SetName(Form("efficiency_%d_%d",iCase,count));
 808   he->SetTitle(Form("efficiency_%d_%d",iCase,count));
 809   return he;
 810 }
 811
 812 THnSparse* CreateResponseMatrix(){
 813   // Creates an empty response matrix as in the Spectrum2 Task
 814
 815   const Int_t kNvar   = 3 ; //number of variables on the grid:pt,eta, phi
 816   const Double_t kPtmin = 0.0, kPtmax = 320.;
 817   const Double_t kEtamin = -3.0, kEtamax = 3.0;
 818   const Double_t kPhimin = 0., kPhimax = 2. * TMath::Pi();
 819   const Double_t kZmin = 0., kZmax = 1;
 820
 821   //arrays for the number of bins in each dimension
 822   Int_t iBin[kNvar];
 823   iBin[0] = 320; //bins in pt
 824   iBin[1] =  1; //bins in eta
 825   iBin[2] = 1; // bins in phi
 826
 827   //arrays for lower bounds :
 828   Double_t* binEdges[kNvar];
 829   for(Int_t ivar = 0; ivar < kNvar; ivar++)
 830     binEdges[ivar] = new Double_t[iBin[ivar] + 1];
 831
 832   //values for bin lower bounds
 833   //  for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)TMath::Power(10,TMath::Log10(kPtmin) + (TMath::Log10(kPtmax)-TMath::Log10(kPtmin))/iBin[0]*(Double_t)i);
 834   for(Int_t i=0; i<=iBin[0]; i++) binEdges[0][i]=(Double_t)kPtmin  + (kPtmax-kPtmin)/(Double_t)iBin[0]*(Double_t)i;
 835   for(Int_t i=0; i<=iBin[1]; i++) binEdges[1][i]=(Double_t)kEtamin  + (kEtamax-kEtamin)/iBin[1]*(Double_t)i;
 836   for(Int_t i=0; i<=iBin[2]; i++) binEdges[2][i]=(Double_t)kPhimin  + (kPhimax-kPhimin)/iBin[2]*(Double_t)i;
 837
 838   Int_t thnDim[2*kNvar];
 839   for (int k=0; k<kNvar; k++) {
 840     //first half  : reconstructed
 841     //second half : MC
 842     thnDim[k]      = iBin[k];
 843     thnDim[k+kNvar] = iBin[k];
 844   }
 845
 846   THnSparseF* fhnCorrelation = new THnSparseF("fhnCorrelation","THnSparse with correlations",2*kNvar,thnDim);
 847   for (int k=0; k<kNvar; k++) {
 848     fhnCorrelation->SetBinEdges(k,binEdges[k]);
 849     fhnCorrelation->SetBinEdges(k+kNvar,binEdges[k]);
 850   }
 851   fhnCorrelation->Sumw2();
 852   return  fhnCorrelation;
 853
 854 }
 855
 856 THnSparse* GetResponseMatrix(TList *fContainer) {
 857   THnSparse* h = (THnSparse*)fContainer->FindObject("fhnCorrelation");
 858   if(fDebug)Printf("Line: %d %s Entries %f",__LINE__,h->GetName(),h->GetEntries());
 859   if(iRebin==1)return h;
 860   Int_t fRebin[6];
 861   fRebin[0] = iRebin;
 862   fRebin[1] = 1;
 863   fRebin[2] = 1;
 864   fRebin[3] = iRebin;
 865   fRebin[4] = 1;
 866   fRebin[5] = 1;
 867
 868   THnSparse* h2 = h->Rebin(fRebin);
 869   return h2;
 870 }
 871
 872 THnSparse* CreateGuessed(const THnSparse* h) {
 873
 874   // is already rebinned
 875   THnSparse* guessed = (THnSparse*) h->Clone();
 876   static int icount = 0;
 877   guessed->SetName(Form("%s_guess_%d",h->GetName(),icount++));
 878
 879
 880   return guessed;
 881
 882 }
 883
 884 void PrintErrors(THnSparse *h){
 885
 886   for(Long_t i = 0; i< h->GetNbins();i++){
 887     Double_t val = h->GetBinContent(i);
 888     Double_t err = h->GetBinError  (i);
 889     if(val>0){
 890       Printf("%s %ld %1.3E +- %1.3E rel Error: %1.3E",h->GetName(),i,val,err,err/val);
 891     }
 892   }
 893
 894 }