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("fHistHadronDepositsAllCent");
 207   TH2F  *fHistHadronDepositsRecoMult = l->FindObject("fHistHadronDepositsRecoCent");
 208   TH2F *eff2D;
 209   if(fHistHadronDepositsRecoMult && fHistHadronDepositsAllMult ){
 210     eff2D = (TH2F*) bayneseffdiv2D(fHistHadronDepositsRecoMult,fHistHadronDepositsAllMult,"eff2D");
 211   }
 212   else{
 213     cerr<<"Warning!  Did not calculate reconstruction efficiency!!"<<endl;
 214     eff2D =  new TH2F("eff2D", "eff2D",200, 0, 10,20,-0.5,19.5);
 215   }
 216   //cor->SetReconstructionEfficiency(eff2D);
 217
 218   TH2F  *fHistGammasGeneratedMult = l->FindObject("fHistGammasGeneratedCent");
 219   TH2F  *fHistGammasFoundMult = l->FindObject("fHistGammasFoundCent");
 220   TH2F *gammaEff2D;
 221   if(fHistGammasGeneratedMult && fHistGammasFoundMult){
 222     gammaEff2D = (TH2F*) bayneseffdiv2D(fHistGammasFoundMult,fHistGammasGeneratedMult,"gammaEff2D");
 223   }
 224   else{
 225     cerr<<"Warning!  Did not calculate reconstruction efficiency!!"<<endl;
 226     gammaEff2D =  new TH2F("gammaEff2D", "gammaEff2D",200, 0, 10,20,-0.5,19.5);
 227   }
 228
 229   AliAnalysisEtTrackMatchCorrections *cor = new AliAnalysisEtTrackMatchCorrections(("TmCorrections"+detector).Data(), fitcharged, fitneutral, fitgamma, fitsecondary, *eff2D,
 230                                                                                    meanCharged, meanNeutral, meanGamma, meanSecondary );
 231
 232
 233   TF1 *func = generateRecEffFunction(p0, p1);
 234   AliAnalysisEtRecEffCorrection *recor = new AliAnalysisEtRecEffCorrection(("ReCorrections"+detector).Data(), *func,*gammaEff2D, 1000);
 235   TFile *outfile = TFile::Open("calocorrections.root","RECREATE");
 236   cor->Write();
 237   recor->Write();
 238   return 0;
 239 }
 240
 241
 242 TF1* generateRecEffFunction(Double_t p0, Double_t p1)
 243 {
 244   TF1 *f = new TF1("receff", "x/([0] + x*[1])", 0, 200);
 245
 246   Double_t params[2] = {p0, p1};
 247   f->SetParameters(params);
 248   return f;
 249 }
 250
 251 TH1* bayneseffdiv(TH1* numerator, TH1* denominator,Char_t* name)
 252 {
 253     if(!numerator){
 254       cerr<<"Error:  numerator does not exist!"<<endl;
 255       return NULL;
 256     }
 257     if(!denominator){
 258       cerr<<"Error:  denominator does not exist!"<<endl;
 259       return NULL;
 260     }
 261     TH1* result = (TH1*)numerator->Clone(name);
 262     Int_t nbins = numerator->GetNbinsX();
 263     for (Int_t ibin=0; ibin<= nbins+1; ++ibin) {
 264       Double_t numeratorVal = numerator->GetBinContent(ibin);
 265       Double_t denominatorVal = denominator->GetBinContent(ibin);
 266       // Check if the errors are right or the thing is scaled
 267       Double_t numeratorValErr = numerator->GetBinError(ibin);
 268       if (!(numeratorValErr==0. || numeratorVal ==0.) ) {
 269         Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal;
 270         numeratorVal /= rescale;
 271       }
 272       Double_t denominatorValErr = denominator->GetBinError(ibin);
 273       if (!(denominatorValErr==0. || denominatorVal==0. )) {
 274         Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal;
 275         denominatorVal /= rescale;
 276       }
 277       Double_t quotient = 0.;
 278       if (denominatorVal!=0.) {
 279         quotient = numeratorVal/denominatorVal;
 280       }
 281       Double_t quotientErr=0;
 282       quotientErr = TMath::Sqrt(
 283                                 (numeratorVal+1.0)/(denominatorVal+2.0)*
 284                                 ((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0)));
 285       result->SetBinContent(ibin,quotient);
 286       result->SetBinError(ibin,quotientErr);
 287       //cout<<"Setting bin "<<ibin<<" to "<<quotient<<" "<<numeratorVal<<"/"<<denominatorVal<<endl;
 288     }
 289     return result;
 290 }
 291
 292
 293 TH2* bayneseffdiv2D(TH2* numerator, TH2* denominator,Char_t* name)
 294 {
 295   if(!numerator){
 296     cerr<<"Error:  numerator does not exist!"<<endl;
 297     return NULL;
 298   }
 299   if(!denominator){
 300     cerr<<"Error:  denominator does not exist!"<<endl;
 301     return NULL;
 302   }
 303   TH2* result = (TH2*)numerator->Clone(name);
 304   Int_t nbinsX = numerator->GetNbinsX();
 305   Int_t nbinsY = numerator->GetNbinsY();
 306   for (Int_t ibin=0; ibin<= nbinsX+1; ++ibin) {
 307     for (Int_t jbin=0; jbin<= nbinsY+1; ++jbin) {
 308       Double_t numeratorVal = numerator->GetBinContent(ibin,jbin);
 309       Double_t denominatorVal = denominator->GetBinContent(ibin,jbin);
 310       // Check if the errors are right or the thing is scaled
 311       Double_t numeratorValErr = numerator->GetBinError(ibin,jbin);
 312       if (!(numeratorValErr==0. || numeratorVal ==0.) ) {
 313         Double_t rescale = numeratorValErr*numeratorValErr/numeratorVal;
 314         numeratorVal /= rescale;
 315       }
 316       Double_t denominatorValErr = denominator->GetBinError(ibin,jbin);
 317       if (!(denominatorValErr==0. || denominatorVal==0. )) {
 318         Double_t rescale = denominatorValErr*denominatorValErr/denominatorVal;
 319         denominatorVal /= rescale;
 320       }
 321       Double_t quotient = 0.;
 322       if (denominatorVal!=0.) {
 323         quotient = numeratorVal/denominatorVal;
 324       }
 325       Double_t quotientErr=0;
 326       quotientErr = TMath::Sqrt(
 327                                 (numeratorVal+1.0)/(denominatorVal+2.0)*
 328                                 ((numeratorVal+2.0)/(denominatorVal+3.0)-(numeratorVal+1.0)/(denominatorVal+2.0)));
 329       result->SetBinContent(ibin,jbin,quotient);
 330       result->SetBinError(ibin,jbin,quotientErr);
 331       //cout<<"Setting bin "<<ibin<<" to "<<quotient<<" "<<numeratorVal<<"/"<<denominatorVal<<endl;
 332     }
 333   }
 334   return result;
 335 }