Modification of AliFemto code to use 4D corrections map
authorlgraczyk <lgraczyk@cern.ch>
Tue, 24 Jun 2014 14:52:01 +0000 (16:52 +0200)
committerlgraczyk <lgraczyk@cern.ch>
Tue, 24 Jun 2014 14:54:20 +0000 (16:54 +0200)
PWGCF/FEMTOSCOPY/AliFemto/AliFemtoEventReaderAOD.cxx
PWGCF/FEMTOSCOPY/AliFemto/AliFemtoTrack.cxx
PWGCF/FEMTOSCOPY/AliFemto/AliFemtoTrack.h
PWGCF/FEMTOSCOPY/AliFemtoUser/AliFemtoCorrFctnDEtaDPhiCorrections.cxx
PWGCF/FEMTOSCOPY/AliFemtoUser/AliFemtoCorrFctnDEtaDPhiCorrections.h

index dd7c235..7fc30bd 100644 (file)
@@ -844,6 +844,7 @@ void AliFemtoEventReaderAOD::CopyAODtoFemtoTrack(AliAODTrack *tAodTrack,
 
   fEvent->GetPrimaryVertex()->GetPosition(fV1);
   //  fEvent->GetPrimaryVertex()->GetXYZ(fV1);
+  tFemtoTrack->SetPrimaryVertex(fV1);
 
   tFemtoTrack->SetCharge(tAodTrack->Charge());
 
index 808a3f0..1e2210d 100644 (file)
@@ -96,6 +96,10 @@ AliFemtoTrack::AliFemtoTrack():
       fNominalTpcPoints[i].SetZ(0);
     }
   //  cout << "Created track " << this << endl;
+
+  fVertex[0] = -9999;
+  fVertex[1] = -9999;
+  fVertex[2] = -9999;
 }
 
 
@@ -244,6 +248,11 @@ AliFemtoTrack::AliFemtoTrack(const AliFemtoTrack& t) :
     //fGlobalEmissionPoint->SetT(t.fGlobalEmissionPoint->e());
     //  cout << "Created track " << this << endl;
   }
+
+
+  fVertex[0] = t.fVertex[0];
+  fVertex[1] = t.fVertex[1];
+  fVertex[2] = t.fVertex[2];
 }
 
 AliFemtoTrack& AliFemtoTrack::operator=(const AliFemtoTrack& aTrack)
@@ -351,6 +360,10 @@ AliFemtoTrack& AliFemtoTrack::operator=(const AliFemtoTrack& aTrack)
     //  cout << "Created track " << this << endl;
   }
 
+  fVertex[0] = aTrack.fVertex[0];
+  fVertex[1] = aTrack.fVertex[1];
+  fVertex[2] = aTrack.fVertex[2];
+
   return *this;
 }
 
@@ -738,3 +751,18 @@ void                   AliFemtoTrack::SetGlobalEmissionPoint(Double_t aRx, Doubl
   }
 }
 //_______________________
+
+
+void AliFemtoTrack::SetPrimaryVertex(double* vertex)
+{
+  fVertex[0] = vertex[0];
+  fVertex[1] = vertex[1];
+  fVertex[2] = vertex[2];
+}
+
+void AliFemtoTrack::GetPrimaryVertex(double* vertex)
+{
+  vertex[0] = fVertex[0];
+  vertex[1] = fVertex[1];
+  vertex[2] = fVertex[2];
+}
index 69f3895..1fb5a6b 100644 (file)
@@ -200,6 +200,10 @@ public:
   void                   SetGlobalEmissionPoint(const AliFemtoThreeVector& aPos);
   void                   SetGlobalEmissionPoint(Double_t aRx, Double_t aRy, Double_t aRz);
  
+
+  void SetPrimaryVertex(double *vertex);
+  void GetPrimaryVertex(double *vertex);
+
   //Alice stuff
   enum {
     kITSin=0x0001,kITSout=0x0002,kITSrefit=0x0004,kITSpid=0x0008,
@@ -288,6 +292,8 @@ public:
   Double_t               fMass;          // True particle mass
   AliFemtoThreeVector   *fGlobalEmissionPoint;
 
+  double fVertex[3];
+
 };
 
 //inline const float* AliFemtoTrack::NSigma() const 
index d422e30..265f693 100644 (file)
@@ -13,6 +13,7 @@
 //#include "AliFemtoHisto.hh"
 #include <cstdio>
 #include <TMath.h>
+#include "THn.h"
 
 #ifdef __ROOT__ 
 ClassImp(AliFemtoCorrFctnDEtaDPhiCorrections)
@@ -52,7 +53,24 @@ AliFemtoCorrFctnDEtaDPhiCorrections::AliFemtoCorrFctnDEtaDPhiCorrections(char* t
   fpTab(0),
   fPartType(kNoCorrection),
   fphiL(0),
-  fphiT(0)
+  fphiT(0),
+  ifileCorrTab(0),
+  fdoPtCorr(0),
+  fdoEtaCorr(0),
+  fdoPhiCorr(0),
+  fdoZVertCorr(0),
+  fpartType1(0),
+  fpartType2(0),
+  fhntReco1(0),
+  fhntReco2(0),
+  fh1Reco1(0),
+  fh1Reco2(0),
+  fh2Reco1(0),
+  fh2Reco2(0),
+  fh3Reco1(0),
+  fh3Reco2(0),
+  fhCont1(0),
+  fhCont2(0)
 {
 
   fphiL = (-(int)(aPhiBins/4)+0.5)*2.*TMath::Pi()/aPhiBins;
@@ -177,7 +195,24 @@ AliFemtoCorrFctnDEtaDPhiCorrections::AliFemtoCorrFctnDEtaDPhiCorrections(const A
   fpTab(0),
   fPartType(kNoCorrection),
   fphiL(0),
-  fphiT(0)
+  fphiT(0),
+  ifileCorrTab(0),
+  fdoPtCorr(0),
+  fdoEtaCorr(0),
+  fdoPhiCorr(0),
+  fdoZVertCorr(0),
+  fpartType1(0),
+  fpartType2(0),
+  fhntReco1(0),
+  fhntReco2(0),
+  fh1Reco1(0),
+  fh1Reco2(0),
+  fh2Reco1(0),
+  fh2Reco2(0),
+  fh3Reco1(0),
+  fh3Reco2(0),
+  fhCont1(0),
+  fhCont2(0)
 {
   // copy constructor
   if (aCorrFctn.fDPhiDEtaNumerator)
@@ -286,6 +321,18 @@ AliFemtoCorrFctnDEtaDPhiCorrections::~AliFemtoCorrFctnDEtaDPhiCorrections(){
   delete fPtCorrectionsDen;
   delete fEtaCorrectionsNum;
   delete fEtaCorrectionsDen;
+
+
+  delete fhntReco1;
+  delete fhntReco2;
+  delete fh1Reco1;
+  delete fh1Reco2;
+  delete fh2Reco1;
+  delete fh2Reco2;
+  delete fh3Reco1;
+  delete fh3Reco2;
+  delete fhCont1;
+  delete fhCont2;
 }
 //_________________________
 AliFemtoCorrFctnDEtaDPhiCorrections& AliFemtoCorrFctnDEtaDPhiCorrections::operator=(const AliFemtoCorrFctnDEtaDPhiCorrections& aCorrFctn)
@@ -437,8 +484,14 @@ void AliFemtoCorrFctnDEtaDPhiCorrections::AddRealPair( AliFemtoPair* pair){
    double pt1 = TMath::Hypot(px1, py1);
    double pt2 = TMath::Hypot(px2, py2);
 
+   double vert1[3];
+   pair->Track1()->Track()->GetPrimaryVertex(vert1);
+   double vert2[3];
+   pair->Track2()->Track()->GetPrimaryVertex(vert2);
+
    double corrweight;
-   if (fIfCorrection) corrweight = CalculateCorrectionWeight(pt1, pt2);
+   //if (fIfCorrection) corrweight = CalculateCorrectionWeight(pt1, pt2);
+  if (fIfCorrection) corrweight = CalculateCorrectionWeight(pt1, pt2, eta1, eta2, phi1, phi2, vert1[2], vert2[2]);
 /*   double ptmin = pt1>pt2 ? pt2 : pt1;
 
    double cosphi = (px1*px2 + py1*py2 + pz1*pz2)/
@@ -517,8 +570,15 @@ void AliFemtoCorrFctnDEtaDPhiCorrections::AddMixedPair( AliFemtoPair* pair){
 //     sqrt((px1*px1 + py1*py1 + pz1*pz1)*(px2*px2 + py2*py2 + pz2*pz2));
 
 
+  double vert1[3];
+   pair->Track1()->Track()->GetPrimaryVertex(vert1);
+   double vert2[3];
+   pair->Track2()->Track()->GetPrimaryVertex(vert2);
+
    double corrweight=-999;
-   if (fIfCorrection) corrweight = CalculateCorrectionWeight(pt1, pt2);
+   //if (fIfCorrection) corrweight = CalculateCorrectionWeight(pt1, pt2);
+  if (fIfCorrection) corrweight = CalculateCorrectionWeight(pt1, pt2, eta1, eta2, phi1, phi2, vert1[2], vert2[2]);
+  
   
    if(fIfCorrection)
       fDPhiDEtaDenominator->Fill(dphi, deta, corrweight);
@@ -659,52 +719,209 @@ void AliFemtoCorrFctnDEtaDPhiCorrections::SetDoCorrectionsHist(CorrectionType do
   fIfCorrectionHist = doCorr;
 }
 
-void AliFemtoCorrFctnDEtaDPhiCorrections::LoadCorrectionTabFromFile(const char *pTtab, const char *corrTab)
+
+
+void AliFemtoCorrFctnDEtaDPhiCorrections::LoadCorrectionTabFromROOTFile(const char *file, ParticleType partType1, ParticleType partType2, bool doPtCorr, bool doEtaCorr, bool doPhiCorr, bool doZVertCorr)
 {
  
-  double val=-10000;
+  ifileCorrTab = TFile::Open(file);
+  fdoPtCorr = doPtCorr;
+  fdoEtaCorr = doEtaCorr;
+  fdoPhiCorr = doPhiCorr;
+  fdoZVertCorr = doZVertCorr;
+  fpartType1 = partType1;
+  fpartType2 = partType2;
+
+  char type1[10];
+  char type2[10];
+
+
+  if(fpartType1==kPion) strcpy(type1,"Pion");
+  else if(fpartType1==kKaon) strcpy(type1,"Kaon");
+  else if (fpartType1==kProton)strcpy(type1,"Proton");
+  else if (fpartType1==kAll) strcpy(type1,"All");
+  else strcpy(type1,"");
+
+  if(fpartType2==kPion) strcpy(type2,"Pion");
+  else if(fpartType2==kKaon) strcpy(type2,"Kaon");
+  else if (fpartType2==kProton) strcpy(type2,"Proton");
+  else if (fpartType2==kAll) strcpy(type2,"All");
+  else strcpy(type1,"");
+
+
+
+  fhntReco1 = (THnT<float>*)(ifileCorrTab->Get(Form("fCorrectionMapData%s",type1)))->Clone();
+  fhntReco2 = (THnT<float>*)(ifileCorrTab->Get(Form("fCorrectionMapData%s",type2)))->Clone();
+  fhCont1 = (TH1D*)(ifileCorrTab->Get(Form("SecondariesContamination%s",type1)))->Clone();
+  fhCont2 = (TH1D*)(ifileCorrTab->Get(Form("SecondariesContamination%s",type2)))->Clone();
 
-  ifstream ifile1;
-  ifile1.open(pTtab);
-  if(ifile1)
+  double fhntReco1_nbins = fhntReco1->GetNbins();
+  double fhntReco2_nbins = fhntReco2->GetNbins();
+
+  int boolSum = fdoPtCorr+fdoEtaCorr+fdoPhiCorr+fdoZVertCorr;
+  /*if(boolSum == 0)
+    {
+      return 1;
+      }*/
+  if(boolSum == 1)
     {
-      int nrEntries1;
-      ifile1>>nrEntries1;
-      fpTab = new double[nrEntries1];
-      int i=0;
-      while(ifile1>>val)
+
+      if(fdoPtCorr == 1)
        {
-         fpTab[i] = val;
-         i++;
+         fh1Reco1 = (TH1F*)(fhntReco1->Projection(0))->Clone();
+         fh1Reco2 = (TH1F*)(fhntReco2->Projection(0))->Clone();
+         fh1Reco1->Scale(1./fhntReco1_nbins*fh1Reco1->GetNbinsX());
+         fh1Reco2->Scale(1./fhntReco2_nbins*fh1Reco2->GetNbinsX());
+           
        }
+
+      else if(fdoEtaCorr == 1)
+       {
+         fh1Reco1 = (TH1F*)(fhntReco1->Projection(1))->Clone();
+         fh1Reco2 = (TH1F*)(fhntReco2->Projection(1))->Clone();
+         fh1Reco1->Scale(1./fhntReco1_nbins*fh1Reco1->GetNbinsX());
+         fh1Reco2->Scale(1./fhntReco2_nbins*fh1Reco2->GetNbinsX());
+       }
+
+      else if(fdoPhiCorr == 1)
+       {
+         fh1Reco1 = (TH1F*)(fhntReco1->Projection(2))->Clone();
+         fh1Reco2 = (TH1F*)(fhntReco2->Projection(2))->Clone();
+         fh1Reco1->Scale(1./fhntReco1_nbins*fh1Reco1->GetNbinsX());
+         fh1Reco2->Scale(1./fhntReco2_nbins*fh1Reco2->GetNbinsX());
+       }
+
+      else if(fdoZVertCorr == 1)
+       {
+         fh1Reco1 = (TH1F*)(fhntReco1->Projection(3))->Clone();
+         fh1Reco2 = (TH1F*)(fhntReco2->Projection(3))->Clone();
+         fh1Reco1->Scale(1./fhntReco1_nbins*fh1Reco1->GetNbinsX());
+         fh1Reco2->Scale(1./fhntReco2_nbins*fh1Reco2->GetNbinsX());
+       }
+
     }
-  else
+
+  else if(boolSum == 2)
     {
-      cout<<"No pT values file open!"<<endl;
+      if(fdoPtCorr == 1 && fdoEtaCorr == 1)
+       {
+         fh2Reco1 = (TH2F*)(fhntReco1->Projection(0,1))->Clone();
+         fh2Reco2 = (TH2F*)(fhntReco2->Projection(0,1))->Clone();        
+         fh2Reco1->Scale(1./fhntReco1_nbins*fh2Reco1->GetNbinsX()*fh2Reco1->GetNbinsY());
+         fh2Reco2->Scale(1./fhntReco2_nbins*fh2Reco2->GetNbinsX()*fh2Reco2->GetNbinsY());
+       }
+
+      if(fdoPtCorr == 1 && fdoPhiCorr == 1)
+       {
+         fh2Reco1 = (TH2F*)(fhntReco1->Projection(0,2))->Clone();
+         fh2Reco2 = (TH2F*)(fhntReco2->Projection(0,2))->Clone();
+         fh2Reco1->Scale(1./fhntReco1_nbins*fh2Reco1->GetNbinsX()*fh2Reco1->GetNbinsY());
+         fh2Reco2->Scale(1./fhntReco2_nbins*fh2Reco2->GetNbinsX()*fh2Reco2->GetNbinsY());       
+       }
+
+      else if(fdoPtCorr == 1 && fdoZVertCorr == 1)
+       {
+         fh2Reco1 = (TH2F*)(fhntReco1->Projection(0,3))->Clone();
+         fh2Reco2 = (TH2F*)(fhntReco2->Projection(0,3))->Clone();
+         fh2Reco1->Scale(1./fhntReco1_nbins*fh2Reco1->GetNbinsX()*fh2Reco1->GetNbinsY());
+         fh2Reco2->Scale(1./fhntReco2_nbins*fh2Reco2->GetNbinsX()*fh2Reco2->GetNbinsY());
+       }
+      else if(fdoEtaCorr == 1 && fdoPhiCorr == 1)
+       {
+         fh2Reco1 = (TH2F*)(fhntReco1->Projection(1,2))->Clone();
+         fh2Reco2 = (TH2F*)(fhntReco2->Projection(1,2))->Clone();
+         fh2Reco1->Scale(1./fhntReco1_nbins*fh2Reco1->GetNbinsX()*fh2Reco1->GetNbinsY());
+         fh2Reco2->Scale(1./fhntReco2_nbins*fh2Reco2->GetNbinsX()*fh2Reco2->GetNbinsY());
+       }
+      else if(fdoEtaCorr == 1 && fdoZVertCorr == 1)
+       {
+         fh2Reco1 = (TH2F*)(fhntReco1->Projection(1,3))->Clone();
+         fh2Reco2 = (TH2F*)(fhntReco2->Projection(1,3))->Clone();
+         fh2Reco1->Scale(1./fhntReco1_nbins*fh2Reco1->GetNbinsX()*fh2Reco1->GetNbinsY());
+         fh2Reco2->Scale(1./fhntReco2_nbins*fh2Reco2->GetNbinsX()*fh2Reco2->GetNbinsY());
+       }
+      else if(fdoPhiCorr == 1 && fdoZVertCorr == 1)
+       {
+         fh2Reco1 = (TH2F*)(fhntReco1->Projection(2,3))->Clone();
+         fh2Reco2 = (TH2F*)(fhntReco2->Projection(2,3))->Clone();
+         fh2Reco1->Scale(1./fhntReco1_nbins*fh2Reco1->GetNbinsX()*fh2Reco1->GetNbinsY());
+         fh2Reco2->Scale(1./fhntReco2_nbins*fh2Reco2->GetNbinsX()*fh2Reco2->GetNbinsY());
+       }
     }
-  ifile1.close();
-  
 
-  ifstream ifile2;
-  ifile2.open(corrTab);
-  if(ifile2)
+
+  else if(boolSum == 3)
     {
-      int nrEntries2;
-      ifile2>>nrEntries2;
-      fCorrFactorTab = new double[nrEntries2];
-      int i=0;
-      while(ifile2>>val)
+      if(fdoPtCorr == 1 && fdoEtaCorr == 1 && fdoPhiCorr == 1)
        {
-         fCorrFactorTab[i] = val;
-         cout<<"fCorrFactorTab: "<<fCorrFactorTab[i]<<endl;
-         i++;
+         fh3Reco1 = (TH3F*)(fhntReco1->Projection(0,1,2))->Clone();
+         fh3Reco2 = (TH3F*)(fhntReco2->Projection(0,1,2))->Clone(); 
+         fh3Reco1->Scale(1./fhntReco1_nbins*fh3Reco1->GetNbinsX()*fh3Reco1->GetNbinsY()*fh3Reco1->GetNbinsZ());
+         fh3Reco2->Scale(1./fhntReco2_nbins*fh3Reco2->GetNbinsX()*fh3Reco2->GetNbinsY()*fh3Reco2->GetNbinsZ());
+       }
+
+      else if(fdoPtCorr == 1 && fdoEtaCorr == 1 && fdoZVertCorr == 1)
+       {
+         fh3Reco1 = (TH3F*)(fhntReco1->Projection(0,1,3))->Clone();
+         fh3Reco2 = (TH3F*)(fhntReco2->Projection(0,1,3))->Clone(); 
+         fh3Reco1->Scale(1./fhntReco1_nbins*fh3Reco1->GetNbinsX()*fh3Reco1->GetNbinsY()*fh3Reco1->GetNbinsZ());
+         fh3Reco2->Scale(1./fhntReco2_nbins*fh3Reco2->GetNbinsX()*fh3Reco2->GetNbinsY()*fh3Reco2->GetNbinsZ());
+       }
+
+      else if(fdoPtCorr == 1 && fdoPhiCorr == 1 && fdoZVertCorr == 1)
+       {
+         fh3Reco1 = (TH3F*)(fhntReco1->Projection(0,2,3))->Clone();
+         fh3Reco2 = (TH3F*)(fhntReco2->Projection(0,2,3))->Clone(); 
+         fh3Reco1->Scale(1./fhntReco1_nbins*fh3Reco1->GetNbinsX()*fh3Reco1->GetNbinsY()*fh3Reco1->GetNbinsZ());
+         fh3Reco2->Scale(1./fhntReco2_nbins*fh3Reco2->GetNbinsX()*fh3Reco2->GetNbinsY()*fh3Reco2->GetNbinsZ());
+       }
+
+      else if(fdoEtaCorr == 1 && fdoPhiCorr == 1 && fdoZVertCorr == 1)
+       {
+         fh3Reco1 = (TH3F*)(fhntReco1->Projection(1,2,3))->Clone();
+         fh3Reco2 = (TH3F*)(fhntReco2->Projection(1,2,3))->Clone(); 
+         fh3Reco1->Scale(1./fhntReco1_nbins*fh3Reco1->GetNbinsX()*fh3Reco1->GetNbinsY()*fh3Reco1->GetNbinsZ());
+         fh3Reco2->Scale(1./fhntReco2_nbins*fh3Reco2->GetNbinsX()*fh3Reco2->GetNbinsY()*fh3Reco2->GetNbinsZ());
        }
     }
-  else
+
+  /*else if(boolSum == 4)
     {
-      cout<<"No corrections file open!"<<endl;
-    }
-  ifile2.close();
+    }*/
+
+  ifileCorrTab->Close();
+
+}
+
+void AliFemtoCorrFctnDEtaDPhiCorrections::LoadCorrectionTabFromROOTFile1D(const char *file, ParticleType partType1, ParticleType partType2)
+{
+  ifileCorrTab = TFile::Open(file);
+
+  fpartType1 = partType1;
+  fpartType2 = partType2;
+
+
+  char type1[10];
+  char type2[10];
+
+
+  if(fpartType1==kPion) strcpy(type1,"Pion");
+  else if(fpartType1==kKaon) strcpy(type1,"Kaon");
+  else if (fpartType1==kProton)strcpy(type1,"Proton");
+  else if (fpartType1==kAll) strcpy(type1,"All");
+  else strcpy(type1,"");
+
+  if(fpartType2==kPion) strcpy(type2,"Pion");
+  else if(fpartType2==kKaon) strcpy(type2,"Kaon");
+  else if (fpartType2==kProton) strcpy(type2,"Proton");
+  else if (fpartType2==kAll) strcpy(type2,"All");
+  else strcpy(type1,"");
+
+  fhCont1 = (TH1D*)(ifileCorrTab->Get(Form("CorrectionFactorPtEffandCont%s",type1)))->Clone();
+  fhCont2 = (TH1D*)(ifileCorrTab->Get(Form("CorrectionFactorPtEffandCont%s",type2)))->Clone();
+  
+  ifileCorrTab->Close();
+
 }
 
 void AliFemtoCorrFctnDEtaDPhiCorrections::SetCorrectionTab(ParticleType partType)
@@ -752,54 +969,320 @@ void AliFemtoCorrFctnDEtaDPhiCorrections::SetCorrectionTab(ParticleType partType
 
 double AliFemtoCorrFctnDEtaDPhiCorrections::CalculateCorrectionWeight(double pT1, double pT2)
 {
+   double w1=0., w2=0.;
+   if(pT1 > fhCont1->GetXaxis()->GetXmin() && pT1 < fhCont1->GetXaxis()->GetXmax() && pT2 > fhCont2->GetXaxis()->GetXmin() && pT2 < fhCont2->GetXaxis()->GetXmax())
+     {
+       w1 = fhCont1->GetBinContent(fhCont1->FindFixBin(pT1));
+       w2 = fhCont2->GetBinContent(fhCont2->FindFixBin(pT2));
+       
+       return w1*w2;
+     } 
+   else
+     return 0;
+}
+
+
+double AliFemtoCorrFctnDEtaDPhiCorrections::CalculateCorrectionWeight(double pT1, double pT2, double eta1, double eta2, double phi1, double phi2, double zvert1, double zvert2)
+{
   
     double w1=0., w2=0.;
-    if(pT1>0 && pT1<5 && pT2>0 && pT2<5)
-    {
-      if(pT1<pT2)
-       {
-          for (int piter = 0 ; piter<200 ; piter++)
-          { 
-        
-             if(pT1>= fpTab[piter] && pT1< fpTab[piter+1])
-              {
-                w1=fCorrFactorTab[piter];
-              }
-             if(pT2>= fpTab[piter] && pT2< fpTab[piter+1])
-             {
-                w2=fCorrFactorTab[piter];
-                break;
-             }
-          }
-       }
-       else if(pT1>pT2)
-       {
-          for (int piter = 0 ; piter<200 ; piter++)
-          {
-             if(pT2>= fpTab[piter] && pT2< fpTab[piter+1])
-                w2=fCorrFactorTab[piter];
-             if(pT1>= fpTab[piter] && pT1< fpTab[piter+1])
-             {
-                w1=fCorrFactorTab[piter];
-                break;
-             }
-          }
-       }
-       else //pT1==pT2
-       {
-          for (int piter = 0 ; piter<200 ; piter++)
-          {
-             if(pT1>= fpTab[piter] && pT1< fpTab[piter+1])
-             {
-                w1=fCorrFactorTab[piter];
-                w2=fCorrFactorTab[piter];
-                break;
-             }
-          }
-       }
-       return w1*w2;
-    }
+    double eps1=0., eps2=0;
+    double cont1=0., cont2=0; //w=(1-cont)/eps
+
+    if(pT1 > fhCont1->GetXaxis()->GetXmin() && pT1 < fhCont1->GetXaxis()->GetXmax() && pT2 > fhCont2->GetXaxis()->GetXmin() && pT2 < fhCont2->GetXaxis()->GetXmax())
+      {
+       cont1 = fhCont1->GetBinContent(fhCont1->FindFixBin(pT1));
+       cont2 = fhCont1->GetBinContent(fhCont2->FindFixBin(pT2));
+      }
     else
-       return 0;
-   return 0;
+      return 0;
+
+    int boolSum = fdoPtCorr+fdoEtaCorr+fdoPhiCorr+fdoZVertCorr;
+    if(boolSum == 0)
+      {
+       return 1;
+      }
+    else if(boolSum == 1)
+      {
+
+       if(fdoPtCorr == 1)
+         {
+           if(pT1 > fh1Reco1->GetXaxis()->GetXmin() && pT1 < fh1Reco1->GetXaxis()->GetXmax() && pT2 > fh1Reco2->GetXaxis()->GetXmin() && pT2 < fh1Reco2->GetXaxis()->GetXmax())
+             {
+               eps1 = fh1Reco1->GetBinContent(fh1Reco1->FindFixBin(pT1));
+               eps2 = fh1Reco2->GetBinContent(fh1Reco2->FindFixBin(pT2));
+               
+               w1 = (1-cont1)/eps1;
+               w2 = (1-cont2)/eps2;
+
+               return w1*w2;
+             }
+           else
+             return 0;     
+         }
+
+       else if(fdoEtaCorr == 1)
+         {
+           if(eta1 > fh1Reco1->GetXaxis()->GetXmin() && eta1 < fh1Reco1->GetXaxis()->GetXmax() && eta2 > fh1Reco2->GetXaxis()->GetXmin() && eta2 < fh1Reco2->GetXaxis()->GetXmax())
+             {
+               eps1 = fh1Reco1->GetBinContent(fh1Reco1->FindFixBin(eta1));
+               eps2 = fh1Reco2->GetBinContent(fh1Reco2->FindFixBin(eta2));
+         
+               w1 = (1-cont1)/eps1;
+               w2 = (1-cont2)/eps2;
+
+               return w1*w2;
+             }
+           else
+             return 0;
+         }
+
+       else if(fdoPhiCorr == 1)
+         {
+           if(phi1 > fh1Reco1->GetXaxis()->GetXmin() && phi1 < fh1Reco1->GetXaxis()->GetXmax() && phi2 > fh1Reco2->GetXaxis()->GetXmin() && phi2 < fh1Reco2->GetXaxis()->GetXmax())
+             {
+               eps1 = fh1Reco1->GetBinContent(fh1Reco1->FindFixBin(phi1));
+               eps2 = fh1Reco2->GetBinContent(fh1Reco2->FindFixBin(phi2));
+
+               w1 = (1-cont1)/eps1;
+               w2 = (1-cont2)/eps2;
+         
+               return w1*w2;
+             }
+           else
+             return 0;
+
+         }
+
+       else if(fdoZVertCorr == 1)
+         {
+           if(zvert1 > fh1Reco1->GetXaxis()->GetXmin() && zvert1 < fh1Reco1->GetXaxis()->GetXmax() && zvert2 > fh1Reco2->GetXaxis()->GetXmin() && zvert2 < fh1Reco2->GetXaxis()->GetXmax())
+             {
+               eps1 = fh1Reco1->GetBinContent(fh1Reco1->FindFixBin(zvert1));
+               eps2 = fh1Reco2->GetBinContent(fh1Reco2->FindFixBin(zvert2));
+         
+               w1 = (1-cont1)/eps1;
+               w2 = (1-cont2)/eps2;
+
+               return w1*w2;
+             }
+           else
+             return 0;
+         }
+
+      }
+
+    else if(boolSum == 2)
+      {
+       if(fdoPtCorr == 1 && fdoEtaCorr == 1)
+         {
+           if(pT1 > fh2Reco1->GetXaxis()->GetXmin() && pT1 < fh2Reco1->GetXaxis()->GetXmax() && pT2 > fh2Reco2->GetXaxis()->GetXmin() && pT2 < fh2Reco2->GetXaxis()->GetXmax() && eta1 > fh2Reco1->GetYaxis()->GetXmin() && eta1 < fh2Reco1->GetYaxis()->GetXmax() && eta2 > fh2Reco2->GetXaxis()->GetXmin() && eta2 < fh2Reco2->GetXaxis()->GetXmax())
+             {
+               eps1 = fh2Reco1->GetBinContent(fh2Reco1->GetXaxis()->FindFixBin(pT1),fh2Reco1->GetYaxis()->FindFixBin(eta1));
+               eps2 = fh2Reco2->GetBinContent(fh2Reco2->GetXaxis()->FindFixBin(pT2),fh2Reco2->GetYaxis()->FindFixBin(eta2));
+
+               w1 = (1-cont1)/eps1;
+               w2 = (1-cont2)/eps2;
+
+               return w1*w2;
+             }
+           else
+             return 0;
+
+         }
+
+       if(fdoPtCorr == 1 && fdoPhiCorr == 1)
+         {
+        
+           if(pT1 > fh2Reco1->GetXaxis()->GetXmin() && pT1 < fh2Reco1->GetXaxis()->GetXmax() && pT2 > fh2Reco2->GetXaxis()->GetXmin() && pT2 < fh2Reco2->GetXaxis()->GetXmax() && phi1 > fh2Reco1->GetYaxis()->GetXmin() && phi1 < fh2Reco1->GetYaxis()->GetXmax() && phi2 > fh2Reco2->GetXaxis()->GetXmin() && phi2 < fh2Reco2->GetXaxis()->GetXmax())
+             {
+               eps1 = fh2Reco1->GetBinContent(fh2Reco1->GetXaxis()->FindFixBin(pT1),fh2Reco1->GetYaxis()->FindFixBin(phi1));
+               eps2 = fh2Reco2->GetBinContent(fh2Reco2->GetXaxis()->FindFixBin(pT2),fh2Reco2->GetYaxis()->FindFixBin(phi2));
+         
+               w1 = (1-cont1)/eps1;
+               w2 = (1-cont2)/eps2;
+
+               return w1*w2;
+             }
+           else
+             return 0;
+
+         }
+
+       else if(fdoPtCorr == 1 && fdoZVertCorr == 1)
+         {
+
+           if(pT1 > fh2Reco1->GetXaxis()->GetXmin() && pT1 < fh2Reco1->GetXaxis()->GetXmax() && pT2 > fh2Reco2->GetXaxis()->GetXmin() && pT2 < fh2Reco2->GetXaxis()->GetXmax() && zvert1 > fh2Reco1->GetYaxis()->GetXmin() && zvert1 < fh2Reco1->GetYaxis()->GetXmax() && zvert2 > fh2Reco2->GetXaxis()->GetXmin() && zvert2 < fh2Reco2->GetXaxis()->GetXmax())
+             {
+               eps1 = fh2Reco1->GetBinContent(fh2Reco1->GetXaxis()->FindFixBin(pT1),fh2Reco1->GetYaxis()->FindFixBin(zvert1));
+               eps2 = fh2Reco2->GetBinContent(fh2Reco2->GetXaxis()->FindFixBin(pT2),fh2Reco2->GetYaxis()->FindFixBin(zvert2));
+
+               w1 = (1-cont1)/eps1;
+               w2 = (1-cont2)/eps2;
+         
+               return w1*w2;
+             }
+           else
+             return 0;
+         }
+       else if(fdoEtaCorr == 1 && fdoPhiCorr == 1)
+         {
+
+           if(eta1 > fh2Reco1->GetXaxis()->GetXmin() && eta1 < fh2Reco1->GetXaxis()->GetXmax() && eta2 > fh2Reco2->GetXaxis()->GetXmin() && eta2 < fh2Reco2->GetXaxis()->GetXmax() && phi1 > fh2Reco1->GetYaxis()->GetXmin() && phi1 < fh2Reco1->GetYaxis()->GetXmax() && phi2 > fh2Reco2->GetXaxis()->GetXmin() && phi2 < fh2Reco2->GetXaxis()->GetXmax())
+             {
+               eps1 = fh2Reco1->GetBinContent(fh2Reco1->GetXaxis()->FindFixBin(eta1),fh2Reco1->GetYaxis()->FindFixBin(phi1));
+               eps2 = fh2Reco2->GetBinContent(fh2Reco2->GetXaxis()->FindFixBin(eta2),fh2Reco2->GetYaxis()->FindFixBin(phi2));
+         
+               w1 = (1-cont1)/eps1;
+               w2 = (1-cont2)/eps2;
+
+               return w1*w2;
+             }
+           else
+             return 0;
+
+
+         }
+       else if(fdoEtaCorr == 1 && fdoZVertCorr == 1)
+         {
+
+           if(eta1 > fh2Reco1->GetXaxis()->GetXmin() && eta1 < fh2Reco1->GetXaxis()->GetXmax() && eta2 > fh2Reco2->GetXaxis()->GetXmin() && eta2 < fh2Reco2->GetXaxis()->GetXmax() && zvert1 > fh2Reco1->GetYaxis()->GetXmin() && zvert1 < fh2Reco1->GetYaxis()->GetXmax() && zvert2 > fh2Reco2->GetXaxis()->GetXmin() && zvert2 < fh2Reco2->GetXaxis()->GetXmax())
+             {
+               eps1 = fh2Reco1->GetBinContent(fh2Reco1->GetXaxis()->FindFixBin(eta1),fh2Reco1->GetYaxis()->FindFixBin(zvert1));
+               eps1 = fh2Reco2->GetBinContent(fh2Reco2->GetXaxis()->FindFixBin(eta2),fh2Reco2->GetYaxis()->FindFixBin(zvert2));
+
+
+               w1 = (1-cont1)/eps1;
+               w2 = (1-cont2)/eps2;
+         
+               return w1*w2;
+             }
+           else
+             return 0;
+
+         }
+       else if(fdoPhiCorr == 1 && fdoZVertCorr == 1)
+         {
+
+           if(phi1 > fh2Reco1->GetXaxis()->GetXmin() && phi1 < fh2Reco1->GetXaxis()->GetXmax() && phi2 > fh2Reco2->GetXaxis()->GetXmin() && phi2 < fh2Reco2->GetXaxis()->GetXmax() && zvert1 > fh2Reco1->GetYaxis()->GetXmin() && zvert1 < fh2Reco1->GetYaxis()->GetXmax() && zvert2 > fh2Reco2->GetYaxis()->GetXmin() && zvert2 < fh2Reco2->GetYaxis()->GetXmax())
+             {
+               eps1 = fh2Reco1->GetBinContent(fh2Reco1->GetXaxis()->FindFixBin(phi1),fh2Reco1->GetYaxis()->FindFixBin(zvert1));
+               eps2 = fh2Reco2->GetBinContent(fh2Reco2->GetXaxis()->FindFixBin(phi2),fh2Reco2->GetYaxis()->FindFixBin(zvert2));
+
+               w1 = (1-cont1)/eps1;
+               w2 = (1-cont2)/eps2;
+         
+               return w1*w2;
+             }
+           else
+             return 0;
+
+         }
+      }
+
+
+    else if(boolSum == 3)
+      {
+       if(fdoPtCorr == 1 && fdoEtaCorr == 1 && fdoPhiCorr == 1)
+         {
+           if(pT1 >fh3Reco1->GetXaxis()->GetXmin() && pT1 <fh3Reco1->GetXaxis()->GetXmax() && 
+               pT2 > fh3Reco2->GetXaxis()->GetXmin() && pT2 <fh3Reco2->GetXaxis()->GetXmax() && 
+               eta1 > fh3Reco1->GetYaxis()->GetXmin() && eta1 <fh3Reco1->GetYaxis()->GetXmax() &&  
+               eta2 > fh3Reco2->GetYaxis()->GetXmin() && eta2 <fh3Reco2->GetYaxis()->GetXmax() &&  
+               phi1 > fh3Reco1->GetZaxis()->GetXmin() && phi1 < fh3Reco1->GetZaxis()->GetXmax() && 
+               phi2 > fh3Reco2->GetZaxis()->GetXmin() && phi2 < fh3Reco2->GetZaxis()->GetXmax())
+             {
+               eps1 = fh3Reco1->GetBinContent(fh3Reco1->GetXaxis()->FindFixBin(pT1),fh3Reco1->GetYaxis()->FindFixBin(eta1),fh3Reco1->GetZaxis()->FindFixBin(phi1));
+               eps2 = fh3Reco2->GetBinContent(fh3Reco2->GetXaxis()->FindFixBin(pT2),fh3Reco2->GetYaxis()->FindFixBin(eta2),fh3Reco2->GetZaxis()->FindFixBin(phi2));
+         
+               w1 = (1-cont1)/eps1;
+               w2 = (1-cont2)/eps2;
+
+               return w1*w2;
+             }
+           else
+             return 0;
+
+
+       }
+
+       else if(fdoPtCorr == 1 && fdoEtaCorr == 1 && fdoZVertCorr == 1)
+         {
+
+           if(pT1 >fh3Reco1->GetXaxis()->GetXmin() && pT1 <fh3Reco1->GetXaxis()->GetXmax() && pT2 > fh3Reco2->GetXaxis()->GetXmin() && pT2 <fh3Reco2->GetXaxis()->GetXmax() && eta1 > fh3Reco1->GetYaxis()->GetXmin() && eta1 <fh3Reco1->GetYaxis()->GetXmax() &&  eta2 > fh3Reco2->GetYaxis()->GetXmin() && eta2 <fh3Reco2->GetYaxis()->GetXmax() &&  zvert1 > fh3Reco1->GetZaxis()->GetXmin() && zvert1 < fh3Reco1->GetZaxis()->GetXmax() &&  zvert2 > fh3Reco2->GetZaxis()->GetXmin() && zvert2 < fh3Reco2->GetZaxis()->GetXmax())
+             {
+               eps1 = fh3Reco1->GetBinContent(fh3Reco1->GetXaxis()->FindFixBin(pT1),fh3Reco1->GetYaxis()->FindFixBin(eta1),fh3Reco1->GetZaxis()->FindFixBin(zvert1));
+               eps2 = fh3Reco2->GetBinContent(fh3Reco2->GetXaxis()->FindFixBin(pT2),fh3Reco2->GetYaxis()->FindFixBin(eta2),fh3Reco2->GetZaxis()->FindFixBin(zvert2));
+         
+               w1 = (1-cont1)/eps1;
+               w2 = (1-cont2)/eps2;
+
+               return w1*w2;
+             }
+           else
+             return 0;
+         }
+
+       else if(fdoPtCorr == 1 && fdoPhiCorr == 1 && fdoZVertCorr == 1)
+         {
+
+           if(pT1 >fh3Reco1->GetXaxis()->GetXmin() && pT1 <fh3Reco1->GetXaxis()->GetXmax() && pT2 > fh3Reco2->GetXaxis()->GetXmin() && pT2 <fh3Reco2->GetXaxis()->GetXmax() && phi1 > fh3Reco1->GetYaxis()->GetXmin() && phi1 <fh3Reco1->GetYaxis()->GetXmax() &&  phi2 > fh3Reco2->GetYaxis()->GetXmin() && phi2 <fh3Reco2->GetYaxis()->GetXmax() &&  zvert1 > fh3Reco1->GetZaxis()->GetXmin() && zvert1 < fh3Reco1->GetZaxis()->GetXmax() &&  zvert2 > fh3Reco2->GetZaxis()->GetXmin() && zvert2 < fh3Reco2->GetZaxis()->GetXmax())
+             {
+               eps1 = fh3Reco1->GetBinContent(fh3Reco1->GetXaxis()->FindFixBin(pT1),fh3Reco1->GetYaxis()->FindFixBin(phi1),fh3Reco1->GetZaxis()->FindFixBin(zvert1));
+               eps2 = fh3Reco2->GetBinContent(fh3Reco2->GetXaxis()->FindFixBin(pT2),fh3Reco2->GetYaxis()->FindFixBin(phi2),fh3Reco2->GetZaxis()->FindFixBin(zvert2));
+         
+               w1 = (1-cont1)/eps1;
+               w2 = (1-cont2)/eps2;
+
+               return w1*w2;
+             }
+           else
+             return 0;
+
+         }
+
+       else if(fdoEtaCorr == 1 && fdoPhiCorr == 1 && fdoZVertCorr == 1)
+         {
+
+           if(eta1 >fh3Reco1->GetXaxis()->GetXmin() && eta1 <fh3Reco1->GetXaxis()->GetXmax() && eta2 > fh3Reco2->GetXaxis()->GetXmin() && eta2 <fh3Reco2->GetXaxis()->GetXmax() && phi1 > fh3Reco1->GetYaxis()->GetXmin() && phi1 <fh3Reco1->GetYaxis()->GetXmax() &&  phi2 > fh3Reco2->GetYaxis()->GetXmin() && phi2 <fh3Reco2->GetYaxis()->GetXmax() &&  zvert1 > fh3Reco1->GetZaxis()->GetXmin() && zvert1 < fh3Reco1->GetZaxis()->GetXmax() &&  zvert2 > fh3Reco2->GetZaxis()->GetXmin() && zvert2 < fh3Reco2->GetZaxis()->GetXmax())
+             {
+               eps1 = fh3Reco1->GetBinContent(fh3Reco1->GetXaxis()->FindFixBin(eta1),fh3Reco1->GetYaxis()->FindFixBin(phi1),fh3Reco1->GetZaxis()->FindFixBin(zvert1));
+               eps2 = fh3Reco2->GetBinContent(fh3Reco2->GetXaxis()->FindFixBin(eta2),fh3Reco2->GetYaxis()->FindFixBin(phi2),fh3Reco2->GetZaxis()->FindFixBin(zvert2));
+         
+               w1 = (1-cont1)/eps1;
+               w2 = (1-cont2)/eps2;
+
+               return w1*w2;
+             }
+           else
+             return 0;
+
+         }
+      }
+
+    else if(boolSum == 4)
+      {
+                                                         
+       if(pT1 > fhntReco1->GetAxis(0)->GetXmin() && pT1 < fhntReco1->GetAxis(0)->GetXmax() && pT2 > fhntReco2->GetAxis(0)->GetXmin() && pT2 < fhntReco2->GetAxis(0)->GetXmax() && eta1 > fhntReco1->GetAxis(1)->GetXmin() && eta1 <fhntReco1->GetAxis(1)->GetXmax() && eta2 > fhntReco2->GetAxis(1)->GetXmin() && eta2 < fhntReco2->GetAxis(1)->GetXmax() && phi1 > fhntReco1->GetAxis(2)->GetXmin() && phi2 < fhntReco2->GetAxis(2)->GetXmax() && phi2 > fhntReco2->GetAxis(2)->GetXmin() && phi2 < fhntReco2->GetAxis(2)->GetXmax() && zvert1 > fhntReco1->GetAxis(3)->GetXmin() && zvert1 < fhntReco1->GetAxis(3)->GetXmax() && zvert2 > fhntReco2->GetAxis(3)->GetXmin() && zvert2 < fhntReco2->GetAxis(3)->GetXmax())
+             {
+
+               int tab1[] = {fhntReco1->GetAxis(0)->FindFixBin(pT1),fhntReco1->GetAxis(1)->FindFixBin(eta1),fhntReco1->GetAxis(2)->FindFixBin(phi1),fhntReco1->GetAxis(3)->FindFixBin(zvert1)};
+               int tab2[] = {fhntReco2->GetAxis(0)->FindFixBin(pT2),fhntReco2->GetAxis(1)->FindFixBin(eta2),fhntReco2->GetAxis(2)->FindFixBin(phi2),fhntReco2->GetAxis(3)->FindFixBin(zvert2)};
+       
+               eps1 = fhntReco1->GetBinContent(tab1);
+               eps2 = fhntReco2->GetBinContent(tab2);
+
+               w1 = (1-cont1)/eps1;
+               w2 = (1-cont2)/eps2;
+                   
+               return w1*w2;
+                 
+             }
+           else
+             return 0;
+
+      }
+    
+    return 0;
+      
 }
index 9d13dcb..b8cdd73 100644 (file)
 
 #include "TH1D.h"
 #include "TH2D.h"
+#include "TH1F.h"
+#include "TH2F.h"
+#include "TH3F.h"
+#include "THn.h"
 #include "THnSparse.h"
 #include "TFile.h"
 #include "AliFemtoCorrFctn.h"
@@ -38,6 +42,9 @@ public:
   void SetDoCorrections(bool doCorr);
   void SetDoCorrectionsHist(CorrectionType doCorr);
   double CalculateCorrectionWeight(double pT1, double pT2);
+  double CalculateCorrectionWeight(double pT1, double pT2, double eta1, double eta2, double phi1, double phi2, double zvert1, double zvert2);
+  void LoadCorrectionTabFromROOTFile1D(const char *file, ParticleType partType1, ParticleType partType2);
+  void LoadCorrectionTabFromROOTFile(const char *file, ParticleType partType1, ParticleType partType2, bool doPtCorr, bool doEtaCorr, bool doPhiCorr, bool doZVertCorr);
   void LoadCorrectionTabFromFile(const char *pTtab, const char *corrTab);
   void SetCorrectionTab(ParticleType partType);
   
@@ -83,10 +90,24 @@ private:
   double fphiL;
   double fphiT;
 
-
-
-
-
+  TFile *ifileCorrTab;
+  bool fdoPtCorr;
+  bool fdoEtaCorr;
+  bool fdoPhiCorr;
+  bool fdoZVertCorr;
+  int fpartType1;
+  int fpartType2;
+
+  THnT<float>* fhntReco1;
+  THnT<float>* fhntReco2;
+  TH1F *fh1Reco1;
+  TH1F *fh1Reco2;
+  TH2F *fh2Reco1;
+  TH2F *fh2Reco2;
+  TH3F *fh3Reco1;
+  TH3F *fh3Reco2;
+  TH1D *fhCont1;
+  TH1D *fhCont2;
 
 #ifdef __ROOT__
   ClassDef(AliFemtoCorrFctnDEtaDPhiCorrections, 1)