]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - PWGLF/totEt/AliAnalysisHadEtCorrections.cxx
- modified addtask for pp for GammaConvV1
[u/mrichter/AliRoot.git] / PWGLF / totEt / AliAnalysisHadEtCorrections.cxx
index 19d0b0b225ac6f64ece019e6d34ed357dd74a4cb..2955ff14c5ae4351c5bcdb4d02f8f0986af5a579 100644 (file)
@@ -73,6 +73,7 @@ AliAnalysisHadEtCorrections::AliAnalysisHadEtCorrections() : TNamed(),
                                                           ,fDataSet(2009)
                                                           ,fProduction("ProductionName")
                                                           ,fProductionDescription("Long production description")
+                                                          ,fSpectraCalcErrorCorrelation(1)
 {//default constructor
   Init();
 }
@@ -168,6 +169,7 @@ AliAnalysisHadEtCorrections::AliAnalysisHadEtCorrections(const AliAnalysisHadEtC
                                                                                              ,fDataSet(g->fDataSet)
                                                                                              ,fProduction(g->fProduction)
                                                                                              ,fProductionDescription(g->fProductionDescription)
+                                                                                             ,fSpectraCalcErrorCorrelation(g->fSpectraCalcErrorCorrelation)
 {//copy constructor
   //SetName(g->GetName());
   fnotIDTPC = new TH1D(*(g->fnotIDTPC));
@@ -200,6 +202,11 @@ Float_t AliAnalysisHadEtCorrections::GetConstantCorrections(Bool_t totEt, Float_
   if(type.Contains("EMCAL")) acceptance = fAcceptanceCorrectionEMCAL;
   if(type.Contains("PHOS")) acceptance = fAcceptanceCorrectionPHOS;
 
+  if(totEt) neutral = fNotHadronicCorrection;
+  else{neutral = fNeutralCorrection;}
+  if(ptcut>0.12){ptcorr = fpTcutCorrectionTPC;}
+  else{ptcorr = fpTcutCorrectionITS;}
+
   if(type.Contains("PiKP")){
     if(ptcut>0.12){ptcorr = fpTcutCorrectionTPC;}
     else{ptcorr = fpTcutCorrectionITS;}
@@ -211,6 +218,7 @@ Float_t AliAnalysisHadEtCorrections::GetConstantCorrections(Bool_t totEt, Float_
       if(ptcut>0.12){ptcorr = ffpTcutCorrectionTPCLow;}
       else{ptcorr = ffpTcutCorrectionITSLow;}
     }
+    neutral = 1.0;
   }
 
   if(type.Contains("High")){//high bound
@@ -229,13 +237,10 @@ Float_t AliAnalysisHadEtCorrections::GetConstantCorrections(Bool_t totEt, Float_
     cout<<"Setting correction factor to "<<correction<<endl;
     return correction;
   }
-
-  if(totEt) neutral = fNotHadronicCorrection;
-  else{neutral = fNeutralCorrection;}
-  if(ptcut>0.12){ptcorr = fpTcutCorrectionTPC;}
-  else{ptcorr = fpTcutCorrectionITS;}
-
+  
+  
   correction = acceptance*neutral*ptcorr;
+  cout<<"correction "<<correction<<" = "<<acceptance<<"*"<<neutral<<"*"<<ptcorr<<endl;
   cout<<"Setting correction factor for ";
   if(totEt) cout<<"total et";
   else{cout<<"hadronic et";}
@@ -243,35 +248,100 @@ Float_t AliAnalysisHadEtCorrections::GetConstantCorrections(Bool_t totEt, Float_
   //cout<<" Acceptance "<<acceptance<<" neutral "<<neutral<<" ptcorr "<<ptcorr<<endl;
   return correction;
 
+}
+void AliAnalysisHadEtCorrections::GetTotalEt(Float_t hadEt,Float_t hadEtErr, Bool_t isTPC, Float_t rawEmEt, Float_t rawEmEtError, Float_t scale, Float_t energyScaleError, Float_t minEt, Float_t minEtError, Float_t nonLinError, Float_t neutronCorr, Float_t neutronError, Float_t hadronCorr, Float_t hadronError, Float_t kaonCorr, Float_t kaonError, Float_t secondaryCorr, Float_t secondaryError,Float_t &totEt, Float_t &totEtError, Float_t &totEtStatError){
+  Float_t neutralFracErr = (fNeutralCorrection - fNeutralCorrectionLow)/fNeutralCorrection;
+  Float_t ptcutFracErr = (ffpTcutCorrectionTPCHigh-fpTcutCorrectionTPC)/fpTcutCorrectionTPC;
+  if(!isTPC) ptcutFracErr = (ffpTcutCorrectionITSHigh-fpTcutCorrectionITS)/fpTcutCorrectionITS;
+  Float_t pidFracErr = (fNotIDConstTPCHigh-fNotIDConstTPC)/fNotIDConstTPC;
+  if(!isTPC) pidFracErr = (fNotIDConstITSHigh-fNotIDConstITS)/fNotIDConstITS;
+  Float_t efficiencyFracErr = (fEfficiencyErrorHigh-fEfficiencyErrorLow)/2.0;
+  Float_t backgroundFracErr = (fBackgroundErrorHigh-fBackgroundErrorLow)/2.0;
+  //Float_t fracerr = TMath::Sqrt(neutralFracErr*neutralFracErr+ptcutFracErr*ptcutFracErr+pidFracErr*pidFracErr+efficiencyFracErr*efficiencyFracErr+backgroundFracErr*backgroundFracErr);
+  Float_t emEt = 0;
+  if(minEt>0){
+    emEt = scale*(rawEmEt-hadronCorr-kaonCorr-neutronCorr-secondaryCorr)/minEt;
+    totEtStatError = TMath::Sqrt(hadEtErr*hadEtErr+TMath::Power(scale*rawEmEtError/minEt,2));
+  }
+  totEt = emEt+hadEt;
+  //first add variance due to had et
+  Float_t variance = (neutralFracErr*neutralFracErr+ptcutFracErr*ptcutFracErr+pidFracErr*pidFracErr+backgroundFracErr*backgroundFracErr)*hadEt*hadEt;
+  //then add variance due to em et
+  variance += (energyScaleError*energyScaleError+minEtError/minEt*minEtError/minEt+nonLinError*nonLinError)*emEt*emEt+scale/minEt*(neutronError*neutronError+hadronError*hadronError+kaonError*kaonError+secondaryError*secondaryError);
+  //add efficiency errors separately - these are basically 100% correlated between Em ET and had Et
+  variance +=efficiencyFracErr*efficiencyFracErr*totEt*totEt;
+  //correlated errors
+//   variance +=fSpectraCalcErrorCorrelation*neutralFracErr*hadEt*(pidFracErr*hadEt+backgroundFracErr*hadEt+scale*kaonError/minEt);
+//   variance +=fSpectraCalcErrorCorrelation*pidFracErr*hadEt*(backgroundFracErr*hadEt+scale*kaonError/minEt+neutralFracErr*hadEt) ;
+//   variance +=fSpectraCalcErrorCorrelation*backgroundFracErr*hadEt*(pidFracErr*hadEt+scale*kaonError/minEt+neutralFracErr*hadEt);
+//   variance +=fSpectraCalcErrorCorrelation*scale*kaonError/minEt*(pidFracErr*hadEt+backgroundFracErr*hadEt+neutralFracErr*hadEt);
+  variance +=fSpectraCalcErrorCorrelation*neutralFracErr*hadEt*(backgroundFracErr*hadEt+scale*kaonError/minEt);
+  variance +=fSpectraCalcErrorCorrelation*backgroundFracErr*hadEt*(scale*kaonError/minEt+neutralFracErr*hadEt);
+  variance +=fSpectraCalcErrorCorrelation*scale*kaonError/minEt*(backgroundFracErr*hadEt+neutralFracErr*hadEt);
+  // cout<<"Correlated error neutral "<<neutralFracErr*hadEt<<" pid "<<pidFracErr*hadEt<<" bkgd "<<backgroundFracErr*hadEt<<" kaon "<<scale*kaonError/minEt<<endl;
+  totEtError = TMath::Sqrt(variance);
+  //cout<<"totEt "<<totEt<<endl;
+
+
 }
 Float_t AliAnalysisHadEtCorrections::GetSystematicErrorBound(Float_t et,Bool_t isLowBound, Bool_t isHadronic, Bool_t isTPC) const{
   //we calculate factors for each value and then multiply them to get the overall bounds
   //neutral corrections, pt cut, pid, efficiency, background
-  float neutral = 1.0;
-  float ptcut = 1.0;
-  float pid = 1.0;
-  float efficiency = 1.0;
-  float background = 1.0;
-  if(isLowBound){//is lower bound
-    if(isHadronic) neutral= fNeutralCorrectionLow/fNeutralCorrection;
-    else{neutral = fNotHadronicCorrectionLow/fNotHadronicCorrection;}
-    if(isTPC) ptcut = ffpTcutCorrectionTPCLow/fpTcutCorrectionTPC;
-    else{ptcut = ffpTcutCorrectionITSLow/fpTcutCorrectionITS;}
-    pid = fNotIDConstTPCLow/fNotIDConstTPC;
-    efficiency = fEfficiencyErrorLow;
-    background = fBackgroundErrorLow;
-  }
-  else{//is higher bound
-    if(isHadronic) neutral= fNeutralCorrectionHigh/fNeutralCorrection;
-    else{neutral= fNotHadronicCorrectionHigh/fNotHadronicCorrection;}
-    if(isTPC) ptcut = ffpTcutCorrectionTPCHigh/fpTcutCorrectionTPC;
-    else{ptcut = ffpTcutCorrectionITSHigh/fpTcutCorrectionITS;}
-    pid = fNotIDConstTPCHigh/fNotIDConstTPC;
-    efficiency = fEfficiencyErrorHigh;
-    background = fBackgroundErrorHigh;
-  }
-  //cout<<"neutral "<<neutral<<" ptcut "<<ptcut<<" pid "<<pid<<" efficiency "<<efficiency<<" background "<<background<<" overall "<<neutral*ptcut*pid*efficiency*background<<endl;
-  return neutral*ptcut*pid*efficiency*background*et;
+
+//   Float_t neutral = fNeutralCorrection;
+//   if(!isHadronic) neutral = fNotHadronicCorrection;
+//   Float_t ptcut = fpTcutCorrectionTPC;
+//   if(!isTPC) ptcut = fpTcutCorrectionITS;
+//   Float_t pid = fNotIDConstTPC;
+//   Float_t efficiency = (fEfficiencyErrorHigh+fEfficiencyErrorLow)/2.0;
+//   Float_t background = fBackground;
+
+  Float_t neutralFracErr = (fNeutralCorrection - fNeutralCorrectionLow)/fNeutralCorrection;
+  if(!isHadronic){
+    neutralFracErr = (fNotHadronicCorrection - fNotHadronicCorrectionLow)/fNotHadronicCorrection;
+  }
+  Float_t ptcutFracErr = (ffpTcutCorrectionTPCHigh-fpTcutCorrectionTPC)/fpTcutCorrectionTPC;
+  if(!isTPC){
+    ptcutFracErr = (ffpTcutCorrectionTPCHigh-fpTcutCorrectionTPC)/fpTcutCorrectionTPC;
+  }
+  Float_t pidFracErr = (fNotIDConstTPCHigh-fNotIDConstTPC)/fNotIDConstTPC;
+  Float_t efficiencyFracErr = (fEfficiencyErrorHigh-fEfficiencyErrorLow)/2.0;
+  Float_t backgroundFracErr = (fBackgroundErrorHigh-fBackgroundErrorLow)/2.0;
+  //Float_t fracerr = TMath::Sqrt(neutralFracErr*neutralFracErr+ptcutFracErr*ptcutFracErr+pidFracErr*pidFracErr+efficiencyFracErr*efficiencyFracErr+backgroundFracErr*backgroundFracErr+fSpectraCalcErrorCorrelation*(neutralFracErr*pidFracErr+neutralFracErr*backgroundFracErr+backgroundFracErr*pidFracErr));
+  Float_t fracerr = TMath::Sqrt(neutralFracErr*neutralFracErr+ptcutFracErr*ptcutFracErr+pidFracErr*pidFracErr+efficiencyFracErr*efficiencyFracErr+backgroundFracErr*backgroundFracErr+fSpectraCalcErrorCorrelation*(neutralFracErr*backgroundFracErr));
+  //cout<<"fracerrs neutral "<<neutralFracErr<<" ptcut "<<ptcutFracErr<<" pid "<<pidFracErr<<" efficiency "<<efficiencyFracErr<<" bkgd "<<backgroundFracErr<<" total fracerr "<<fracerr<<endl;
+  if(isLowBound){
+    return et*(1.0-fracerr);
+  }
+  else{
+    return et*(1.0+fracerr);
+  }
+
+//   Float_t neutral = 1.0;
+//   Float_t ptcut = 1.0;
+//   Float_t pid = 1.0;
+//   Float_t efficiency = 1.0;
+//   Float_t background = 1.0;
+//   if(isLowBound){//is lower bound
+//     if(isHadronic) neutral= fNeutralCorrectionLow/fNeutralCorrection;
+//     else{neutral = fNotHadronicCorrectionLow/fNotHadronicCorrection;}
+//     if(isTPC) ptcut = ffpTcutCorrectionTPCLow/fpTcutCorrectionTPC;
+//     else{ptcut = ffpTcutCorrectionITSLow/fpTcutCorrectionITS;}
+//     pid = fNotIDConstTPCLow/fNotIDConstTPC;
+//     efficiency = fEfficiencyErrorLow;
+//     background = fBackgroundErrorLow;
+//   }
+//   else{//is higher bound
+//     if(isHadronic) neutral= fNeutralCorrectionHigh/fNeutralCorrection;
+//     else{neutral= fNotHadronicCorrectionHigh/fNotHadronicCorrection;}
+//     if(isTPC) ptcut = ffpTcutCorrectionTPCHigh/fpTcutCorrectionTPC;
+//     else{ptcut = ffpTcutCorrectionITSHigh/fpTcutCorrectionITS;}
+//     pid = fNotIDConstTPCHigh/fNotIDConstTPC;
+//     efficiency = fEfficiencyErrorHigh;
+//     background = fBackgroundErrorHigh;
+//   }
+//   //cout<<"neutral "<<neutral<<" ptcut "<<ptcut<<" pid "<<pid<<" efficiency "<<efficiency<<" background "<<background<<" overall "<<neutral*ptcut*pid*efficiency*background<<endl;
+//   return neutral*ptcut*pid*efficiency*background*et;
 }
 // AliAnalysisHadEtCorrections & operator = (const AliAnalysisHadEtCorrections & g) {
 
@@ -611,10 +681,10 @@ void AliAnalysisHadEtCorrections::Report(){//Gives a report on the status of all
     break;
   case 2012:
     cout<<"p+p collisions at 8 TeV"<<endl;
+    break;
   case 2013:
     cout<<"p+Pb collisions at 5 TeV"<<endl;
     break;
-    break;
   case 20111:
     cout<<"p+p collisions at 2.76 TeV"<<endl;
     break;