PWGLF/totEt/macros/calculateCorrections.C

   1 //Oystein Djuvsland
   2 //macro for plotting correction factors for EM et for use with AliAnalysisEt and writes them to a file
   3 #include "TTree.h"
   4 #include "TFile.h"
   5 #include <TList.h>
   6 #include <Rtypes.h>
   7 #include "TCanvas.h"
   8 #include <TH2I.h>
   9 #include <TH2F.h>
  10 #include <TProfile.h>
  11 #include <TFitResultPtr.h>
  12 #include <TFitResult.h>
  13 #include "TF1.h"
  14 #include <iostream>
  15 #include <AliAnalysisEtTrackMatchCorrections.h>
  16 #include <AliAnalysisEtRecEffCorrection.h>
  17
  18 TCanvas *c1 = 0;
  19
  20 TH2* bayneseffdiv2D(TH2* numerator, TH2* denominator,Char_t* name) ;
  21 TH1* bayneseffdiv(TH1* numerator, TH1* denominator,Char_t* name) ;
  22
  23 //creates an empty set of correction factors for debugging purposes
  24 TF1 * generateRecEffFunction(Double_t p0, Double_t p1);
  25 int createDummy(Double_t p0 = 1.0, Double_t p1 = 0.0, char *det = "Phos")
  26 {
  27   TFile *outfile = TFile::Open("calocorrections.root","RECREATE");
  28   TF1 fitneutral("fitneutral","0", 0, 100);
  29   TF1 fitcharged("fitcharged","0", 0, 100);
  30   TF1 fitgamma("fitgamma","0", 0, 100);
  31   TF1 fitsecondary("fitsecondary","0", 0, 100);
  32   AliAnalysisEtTrackMatchCorrections *cor = new AliAnalysisEtTrackMatchCorrections(Form("TmCorrections%s",det), fitcharged, fitneutral, fitgamma, fitsecondary,0,0,0,0);
  33   TF1 *func = generateRecEffFunction(p0, p1);
  34   AliAnalysisEtRecEffCorrection *recor = new AliAnalysisEtRecEffCorrection(Form("ReCorrections%s",det), *func, 1000);
  35
  36   cor->Write();
  37   recor->Write();
  38   outfile->Close();
  39
  40
  41 }
  42 //p0 = correction factors for efficiency from track matching
  43 //p1 = correction factors for efficiency from track matching
  44 //determined from fit of efficiency from track matching, 0.366 from simulations, fit as a function of energy
  45 //p0 is the constant
  46 //p1 is linear term
  47 int calculateCorrections(TString filename="rootFiles/LHC11a10a_bis/Et.ESD.simPbPb.EMCal.LHC11a10a_bis.Run139465.root", Double_t p0 = 0.366, Double_t p1 = 0.0, char *det = "Emcal", Bool_t isSim = kFALSE)
  48 {
  49   TString detector = det;
  50   TString emcal = "Emcal";
  51   c1 = new TCanvas;
  52
  53   TFile *f = TFile::Open(filename, "READ");
  54
  55   TList *l = (TList*)f->Get("out1");
  56
  57   TTree *primaryTree = (TTree*)l->FindObject(("fPrimaryTree"+detector+"MC").Data());
  58   TTree *recTree = (TTree*)l->FindObject(("fEventSummaryTree"+detector+"Rec").Data());
  59   TTree *mcTree = (TTree*)l->FindObject(("fEventSummaryTree"+detector+"MC").Data());
  60   std::cout << primaryTree << " " << recTree << " " << mcTree << std::endl;
  61
  62   Int_t clusterMult = 0;
  63   Int_t nChargedNotRemoved = 0;
  64   Int_t nNeutralNotRemoved = 0;
  65   Int_t nGammaRemoved = 0;
  66   Int_t nSecNotRemoved = 0;
  67
  68   Double_t emEtRec = 0;
  69   Double_t emEtMc = 0;
  70
  71   recTree->SetBranchAddress("fNeutralMultiplicity", &clusterMult);
  72   mcTree->SetBranchAddress("fChargedNotRemoved", &nChargedNotRemoved);
  73   mcTree->SetBranchAddress("fNeutralNotRemoved", &nNeutralNotRemoved);
  74   mcTree->SetBranchAddress("fGammaRemoved", &nGammaRemoved);
  75   mcTree->SetBranchAddress("fSecondaryNotRemoved", &nSecNotRemoved);
  76
  77   Float_t maxMult = 99.5;
  78   Int_t nbins = 100;
  79 //   if(detector==emcal){
  80 //     //100 seems to be sufficient...
  81 //     maxMult = 299.5;
  82 //     nbins = 300;
  83 //  }
  84   TH2I *hChargedVsClusterMult = new TH2I("hChVsMult", "hChVsMult", nbins, -0.5, maxMult, nbins, -0.5, maxMult);
  85   TH2I *hNeutralVsClusterMult = new TH2I("hNeutVsMult", "hNeutVsMult", nbins, -0.5, maxMult, nbins, -0.5, maxMult);
  86   TH2I *hGammaVsClusterMult = new TH2I("hGammaVsMult", "hGammaVsMult", nbins, -0.5, maxMult, nbins, -0.5, maxMult);
  87   TH2I *hSecVsClusterMult = new TH2I("hSecVsMult", "hSecVsMult", nbins, -0.5, maxMult, nbins, -0.5, maxMult);
  88
  89   Int_t nEvents = mcTree->GetEntriesFast();
  90     for(Int_t i = 0; i < nEvents; i++)
  91     {
  92       mcTree->GetEvent(i);
  93       recTree->GetEvent(i);
  94       hChargedVsClusterMult->Fill(clusterMult, nChargedNotRemoved);
  95       hNeutralVsClusterMult->Fill(clusterMult, nNeutralNotRemoved);
  96       hGammaVsClusterMult->Fill(clusterMult, nGammaRemoved);
  97       hSecVsClusterMult->Fill(clusterMult, nSecNotRemoved);
  98     }
  99
 100   c1->Divide(2,2);
 101   c1->cd(1);
 102   TString title = "Charged particles not removed";
 103   TString xtitle = "Cluster Multiplicity";
 104   TString ytitle = "N_{ch}";
 105   hChargedVsClusterMult->SetTitle(title);
 106   hChargedVsClusterMult->GetYaxis()->SetTitle(ytitle);
 107   hChargedVsClusterMult->GetXaxis()->SetTitle(xtitle);
 108   hChargedVsClusterMult->SetStats(0);
 109   hChargedVsClusterMult->Draw();
 110   TProfile *chProf = hChargedVsClusterMult->ProfileX();
 111   chProf->SetStats(0);
 112   chProf->Draw("SAME");
 113   //TF1 fitcharged("fitcharged","([0] + [1]*x)*(0.48/([2] + [3]*[2]))", 0, 100);
 114   TF1 fitcharged("fitcharged","pol2", 0, 100);//fit of number of charged tracks vs detector multiplicity
 115   //if straight line track matching roughly not dependent on centrality
 116 //   fitcharged.FixParameter(2, p0);
 117 //   fitcharged.FixParameter(3, p1);
 118   TFitResultPtr chRes = chProf->Fit(&fitcharged,"S");
 119   TArrayD ch;
 120   if(!chRes)
 121   {
 122     ch = TArrayD(chRes->NPar(),chRes->GetParams());
 123   }
 124   else
 125   {
 126     std::cout << "Could not extract charged contribution params" << std::endl;
 127   }
 128   c1->cd(2);
 129   title = "Neutral particles not removed";
 130   ytitle = "N_{neutral}";
 131   hNeutralVsClusterMult->SetTitle(title);
 132   hNeutralVsClusterMult->GetXaxis()->SetTitle(xtitle);
 133   hNeutralVsClusterMult->GetYaxis()->SetTitle(ytitle);
 134   hNeutralVsClusterMult->SetStats(0);
 135   hNeutralVsClusterMult->Draw();
 136   TProfile *neuProf = hNeutralVsClusterMult->ProfileX();
 137   neuProf->SetStats(0);
 138   neuProf->Draw("SAME");
 139   TF1 fitneutral("fitneutral","pol2", 0, 100);//fit of number of neutral particles in calo that we call background in calo vs detector multiplicity
 140   //may include K0S
 141   TFitResultPtr neuRes = neuProf->Fit(&fitneutral,"S");
 142   TArrayD neu;
 143   if(!neuRes)
 144   {
 145     neu = TArrayD(neuRes->NPar(),neuRes->GetParams());
 146   }
 147   else
 148   {
 149     std::cout << "Could not extract charged contribution params" << std::endl;
 150   }
 151   c1->cd(3);
 152
 153   title = "Gammas removed";
 154   ytitle = "N_{#gamma}";
 155   hGammaVsClusterMult->SetTitle(title);
 156   hGammaVsClusterMult->GetYaxis()->SetTitle(ytitle);
 157   hGammaVsClusterMult->GetXaxis()->SetTitle(xtitle);
 158   hGammaVsClusterMult->SetStats(0);
 159   hGammaVsClusterMult->Draw();
 160   TProfile *gamProf = hGammaVsClusterMult->ProfileX();
 161   gamProf->SetStats(0);
 162   gamProf->Draw("SAME");
 163   TF1 fitgamma("fitgamma","pol2", 0, 100);//fit of number of gammas removed erroneously vs detector multiplicity
 164   TFitResultPtr gammaRes = gamProf->Fit(&fitgamma,"S");
 165   TArrayD gamma;
 166   if(!gammaRes)
 167   {
 168     gamma = TArrayD(gammaRes->NPar(),gammaRes->GetParams());
 169   }
 170   else
 171   {
 172     std::cout << "Could not extract charged contribution params" << std::endl;
 173   }
 174
 175   c1->cd(4);
 176
 177   title = "Secondaries not removed";
 178   ytitle = "N_{sec}";
 179   hSecVsClusterMult->SetTitle(title);
 180   hSecVsClusterMult->GetXaxis()->SetTitle(xtitle);
 181   hSecVsClusterMult->GetYaxis()->SetTitle(ytitle);
 182   hSecVsClusterMult->SetStats(0);
 183   hSecVsClusterMult->Draw();
 184   TProfile *secProf = hSecVsClusterMult->ProfileX();
 185   secProf->SetStats(0);
 186   secProf->Draw("SAME");
 187   TF1 fitsecondary("fitsecondary","pol2", 0, 100);//fit of number of secondary particles that leave deposits in calo vs detector multiplicity
 188   TFitResultPtr secRes = secProf->Fit(&fitsecondary,"S");
 189   TArrayD sec;
 190   if(!secRes)
 191   {
 192     sec = TArrayD(secRes->NPar(),secRes->GetParams());
 193   }
 194   else
 195   {
 196     std::cout << "Could not extract charged contribution params" << std::endl;
 197   }
 198   //ugly hack for getting the energy
 199   //changing number of particles to energy of particles by multiplying by the <energy> and dividing by the efficiency for the average energy
 200   //average energy from each of these:  in excel file
 201   Double_t meanCharged = 0.48/(p0 + p1*0.48);
 202   Double_t meanNeutral = 0.53/(p0 + p1*0.53);
 203   Double_t meanGamma = 0.51/(p0 + p1*0.51);
 204   Double_t meanSecondary = meanGamma;
 205
 206   TH2F  *fHistHadronDepositsAllMult = l->FindObject("fHistHadronDepositsAllMult");
 207   TH2F  *fHistHadronDepositsRecoMult = l->FindObject("fHistHadronDepositsRecoMult");
 208   TH2F *eff2D = (TH2F*) bayneseffdiv2D(fHistHadronDepositsRecoMult,fHistHadronDepositsAllMult,"eff2D");
 209   //cor->SetReconstructionEfficiency(eff2D);
 210
 211   TH2F  *fHistGammasGeneratedMult = l->FindObject("fHistGammasGeneratedMult");
 212   TH2F  *fHistGammasFoundMult = l->FindObject("fHistGammasFoundMult");
 213   TH2F *gammaEff2D = (TH2F*) bayneseffdiv2D(fHistGammasFoundMult,fHistGammasGeneratedMult,"gammaEff2D");
 214
 215   AliAnalysisEtTrackMatchCorrections *cor = new AliAnalysisEtTrackMatchCorrections(("TmCorrections"+detector).Data(), fitcharged, fitneutral, fitgamma, fitsecondary, *eff2D,
 216                                                                                    meanCharged, meanNeutral, meanGamma, meanSecondary );
 217
 218
 219   TF1 *func = generateRecEffFunction(p0, p1);
 220   AliAnalysisEtRecEffCorrection *recor = new AliAnalysisEtRecEffCorrection(("ReCorrections"+detector).Data(), *func,*gammaEff2D, 1000);
 221   TFile *outfile = TFile::Open("calocorrections.root","RECREATE");
 222   cor->Write();
 223   recor->Write();
 224   return 0;
 225 }
 226
 227
 228 TF1* generateRecEffFunction(Double_t p0, Double_t p1)
 229 {
 230   TF1 *f = new TF1("receff", "x/([0] + x*[1])", 0, 200);
 231
 232   Double_t params[2] = {p0, p1};
 233   f->SetParameters(params);
 234   return f;
 235 }
 236
 237 TH1* bayneseffdiv(TH1* numerator, TH1* denominator,Char_t* name)
 238 {
 239     if(!numerator){
 240       cerr<<"Error:  numerator does not exist!"<<endl;
 241       return NULL;
 242     }
 243     if(!denominator){
 244       cerr<<"Error:  denominator does not exist!"<<endl;
 245       return NULL;
 246     }
 247     TH1* result = (TH1*)numerator->Clone(name);
 248     Int_t nbins = numerator->GetNbinsX();
 249     for (Int_t ibin=0; ibin<= nbins+1; ++ibin) {
 250       Double_t numeratorVal = numerator->GetBinContent(ibin);
 251       Double_t denominatorVal = denominator->GetBinContent(ibin);
 252       // Check if the errors are right or the thing is scaled
 253       Double_t numeratorValErr = numerator->GetBinError(ibin);
 254       if (!(numeratorValErr==0. || numeratorVal ==0.) ) {
 255         Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal;
 256         numeratorVal /= rescale;
 257       }
 258       Double_t denominatorValErr = denominator->GetBinError(ibin);
 259       if (!(denominatorValErr==0. || denominatorVal==0. )) {
 260         Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal;
 261         denominatorVal /= rescale;
 262       }
 263       Double_t quotient = 0.;
 264       if (denominatorVal!=0.) {
 265         quotient = numeratorVal/denominatorVal;
 266       }
 267       Double_t quotientErr=0;
 268       quotientErr = TMath::Sqrt(
 269                                 (numeratorVal+1.0)/(denominatorVal+2.0)*
 270                                 ((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0)));
 271       result->SetBinContent(ibin,quotient);
 272       result->SetBinError(ibin,quotientErr);
 273       //cout<<"Setting bin "<<ibin<<" to "<<quotient<<" "<<numeratorVal<<"/"<<denominatorVal<<endl;
 274     }
 275     return result;
 276 }
 277
 278
 279 TH2* bayneseffdiv2D(TH2* numerator, TH2* denominator,Char_t* name)
 280 {
 281   if(!numerator){
 282     cerr<<"Error:  numerator does not exist!"<<endl;
 283     return NULL;
 284   }
 285   if(!denominator){
 286     cerr<<"Error:  denominator does not exist!"<<endl;
 287     return NULL;
 288   }
 289   TH2* result = (TH2*)numerator->Clone(name);
 290   Int_t nbinsX = numerator->GetNbinsX();
 291   Int_t nbinsY = numerator->GetNbinsY();
 292   for (Int_t ibin=0; ibin<= nbinsX+1; ++ibin) {
 293     for (Int_t jbin=0; jbin<= nbinsY+1; ++jbin) {
 294       Double_t numeratorVal = numerator->GetBinContent(ibin,jbin);
 295       Double_t denominatorVal = denominator->GetBinContent(ibin,jbin);
 296       // Check if the errors are right or the thing is scaled
 297       Double_t numeratorValErr = numerator->GetBinError(ibin,jbin);
 298       if (!(numeratorValErr==0. || numeratorVal ==0.) ) {
 299         Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal;
 300         numeratorVal /= rescale;
 301       }
 302       Double_t denominatorValErr = denominator->GetBinError(ibin,jbin);
 303       if (!(denominatorValErr==0. || denominatorVal==0. )) {
 304         Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal;
 305         denominatorVal /= rescale;
 306       }
 307       Double_t quotient = 0.;
 308       if (denominatorVal!=0.) {
 309         quotient = numeratorVal/denominatorVal;
 310       }
 311       Double_t quotientErr=0;
 312       quotientErr = TMath::Sqrt(
 313                                 (numeratorVal+1.0)/(denominatorVal+2.0)*
 314                                 ((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0)));
 315       result->SetBinContent(ibin,jbin,quotient);
 316       result->SetBinError(ibin,jbin,quotientErr);
 317       //cout<<"Setting bin "<<ibin<<" to "<<quotient<<" "<<numeratorVal<<"/"<<denominatorVal<<endl;
 318     }
 319   }
 320   return result;
 321 }