[u/mrichter/AliRoot.git] / PWGLF / FORWARD / analysis2 / AliDisplacedVertexSelection.cxx

#include "AliDisplacedVertexSelection.h"
#include <iostream>
#include <TROOT.h>
#include "AliESDEvent.h"
#include "AliESDZDC.h"
ClassImp(AliDisplacedVertexSelection)
#if 0
; // This is for Emacs
#endif

//____________________________________________________________________
AliDisplacedVertexSelection::AliDisplacedVertexSelection()
  : TObject(), 
    fVertexZ(9999), 
    fCent(100)
{
}
//____________________________________________________________________
AliDisplacedVertexSelection::AliDisplacedVertexSelection(const AliDisplacedVertexSelection& o)
  : TObject(o), 
    fVertexZ(9999), 
    fCent(100)
{
}
//____________________________________________________________________
AliDisplacedVertexSelection&
AliDisplacedVertexSelection::operator=(const AliDisplacedVertexSelection& o)
{
  if (&o == this) return *this; 
  return *this;
}

//____________________________________________________________________
void
AliDisplacedVertexSelection::Output(TList* /*l*/, const char* /* name*/) const
{
}
  
//____________________________________________________________________
void
AliDisplacedVertexSelection::Print(Option_t*) const
{
#if 0
  char ind[gROOT->GetDirLevel()+1];
  for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
  ind[gROOT->GetDirLevel()] = '\0';
  std::cout << std::boolalpha 
	    << std::noboolalpha << std::endl;
#endif
}

//____________________________________________________________________
Bool_t
AliDisplacedVertexSelection::Process(const AliESDEvent* esd)
{
  fVertexZ = 9999; // Default vertex value 
  fCent    = 100;  // Default centrality value 

  // Some constants 
  const Float_t kZDCrefSum        = -66.5;
  const Float_t kZDCrefDelta      =  -2.10;
  const Float_t kZDCsigmaSum      =   3.25;
  const Float_t kZDCsigmaDelta    =   2.25;
  const Float_t kZDCsigmaSumSat   =   2.50;
  const Float_t kZDCsigmaDeltaSat =   1.35;  
  
  // Corrections for magnetic fields
  const Float_t kZEMcorrPlusPlus[21]   = {0.6840,0.7879,0.8722,
					  0.9370,0.9837,1.0137,
					  1.0292,1.0327,1.0271,
					  1.0152,1.0000,0.9844,
					  0.9714,0.9634,0.9626,
					  0.9708,0.9891,1.0181,
					  1.0574,1.1060,1.1617};
  const Float_t kZEMcorrMoinsMoins[21] = {0.7082,0.8012,0.8809,
					  0.9447,0.9916,1.0220,
					  1.0372,1.0395,1.0318,
					  1.0174,1.0000,0.9830,
					  0.9699,0.9635,0.9662,
					  0.9794,1.0034,1.0371,
					  1.0782,1.1224,1.1634};

  // --- Get the ESD object ------------------------------------------
  AliESDZDC* esdZDC = esd->GetESDZDC();
  if (!esdZDC) { 
    AliWarning("No ZDC ESD object!");
    return false; 
  }
  Double_t currentL3   = esd->GetCurrentL3();
  Double_t currentDipo = esd->GetCurrentDip();

  // --- ZDC and ZEM energy signal -----------------------------------
  Double_t zdcEn      = (esdZDC->GetZDCN1Energy()+
			 esdZDC->GetZDCP1Energy()+
			 esdZDC->GetZDCN2Energy()+
			 esdZDC->GetZDCP2Energy());
  Double_t zemEn      = (esdZDC->GetZDCEMEnergy(0)+
			 esdZDC->GetZDCEMEnergy(1))/8.;

  // --- Calculate DeltaT and sumT -----------------------------------
  Double_t deltaTdc = 0;
  Double_t sumTdc   = 0;
   
  for(Int_t i = 0; i < 4; ++i) {
    if(esdZDC->GetZDCTDCData(10,i) != 0) {
      Double_t  tdcCnoCorr = 0.025*(esdZDC->GetZDCTDCData(10,i)-
				    esdZDC->GetZDCTDCData(14,i)); 
      Double_t  tdcC       = esdZDC->GetZDCTDCCorrected(10,i); 
      for(Int_t j = 0; j < 4; ++j) {
	if(esdZDC->GetZDCTDCData(12,j) != 0) {
	  Double_t   tdcAnoCorr = 0.025*(esdZDC->GetZDCTDCData(12,j)-
					 esdZDC->GetZDCTDCData(14,j));
	  Double_t   tdcA       = esdZDC->GetZDCTDCCorrected(12,j);
	  if(esdZDC->TestBit(AliESDZDC::kCorrectedTDCFilled)) {	    
	    deltaTdc = tdcC-tdcA;
	    sumTdc   = tdcC+tdcA;
	  }
	  else {
	    deltaTdc = tdcCnoCorr-tdcAnoCorr;
	    sumTdc   = tdcCnoCorr+tdcAnoCorr;
	  }
	}
      }
    }
  }

  // --- Find the vertex ---------------------------------------------
  if(deltaTdc!=0. || sumTdc!=0.) {
    for (Int_t k = -10; k <= 10; ++k) {
      Float_t zsat  = 2.55F * k;
      Float_t delta = (k == 0 ? kZDCsigmaDelta : kZDCsigmaDeltaSat);
      Float_t sum   = (k == 0 ? kZDCsigmaSum   : kZDCsigmaSumSat);
      Float_t dT    = deltaTdc - kZDCrefDelta - zsat;
      Float_t sT    = sumTdc  - kZDCrefSum - zsat;
      Float_t check = dT * dT / delta / delta + sT * sT / sum  / sum;
      if (check > 1.0 || k == 0) continue; 
      
      // Set the vertex 
      fVertexZ = 37.5 * k;
      
      // Correct zem energy 
      if(currentDipo>0 && currentL3>0) zemEn /= kZEMcorrPlusPlus[k+10];
      if(currentDipo<0 && currentL3<0) zemEn /= kZEMcorrMoinsMoins[k+10];
    }
  }

  // --- Calculate the centrality ------------------------------------
  Float_t c1, c2, c3, c4, c5;
  Int_t runnumber = esd->GetRunNumber();
  if (runnumber < 137722 && runnumber >= 137161 ) {
    c1 = 16992;
    c2 = 345;
    c3 = 2.23902e+02;
    c4 = 1.56667;
    c5 = 9.49434e-05;
  }
  else if (runnumber < 139172 && runnumber >= 137722) {
    c1  = 15000;
    c2  = 295;
    c3  = 2.23660e+02;
    c4 = 1.56664;
    c5 = 8.99571e-05;
  }
  else if (runnumber >= 139172) {
    c1 = 16992;
    c2 = 345;
    c3 = 2.04789e+02;
    c4 = 1.56629;
    c5 = 1.02768e-04;
  }
  else {
    c1 = 16992;
    c2 = 345;
    c3 = 2.04789e+02;
    c4 = 1.56629;
    c5 = 1.02768e-04;
  }

  if (zemEn > c2) { 
    Float_t slope   = (zdcEn + c1) / (zemEn - c2) + c3;
    Float_t zdcCent = (TMath::ATan(slope) - c4) / c5;
    if (zdcCent >= 0) fCent = zdcCent;
  }

  return true;
}

//____________________________________________________________________
Double_t
AliDisplacedVertexSelection::CheckDisplacedVertex(const AliESDEvent* esd) const {
  // Code taken directly from M.Guilbaud.
  AliWarning("This method is deprecated");
  Double_t zvtx = 9999.;
  
  Float_t kZDCrefSum        =-66.5;
  Float_t kZDCrefDelta      =-2.10;
  Float_t kZDCsigmaSum      = 3.25;
  Float_t kZDCsigmaDelta    = 2.25;
  Float_t kZDCsigmaSumSat   = 2.50;
  Float_t kZDCsigmaDeltaSat = 1.35;  
  
  AliESDZDC* esdZDC = esd -> GetESDZDC();
  
  /* Double_t zdcNCEnergy  = esdZDC->GetZDCN1Energy();
     Double_t fZpcEnergy  = esdZDC->GetZDCP1Energy();
     Double_t fZnaEnergy  = esdZDC->GetZDCN2Energy();
     Double_t fZpaEnergy  = esdZDC->GetZDCP2Energy();
     Double_t fZem1Energy = esdZDC->GetZDCEMEnergy(0);
     Double_t fZem2Energy = esdZDC->GetZDCEMEnergy(1);*/

  //Double_t fzdcEn      = (esdZDC->GetZDCN1Energy()+esdZDC->GetZDCP1Energy()+esdZDC->GetZDCN2Energy()+esdZDC->GetZDCP2Energy());
  //Double_t fzemEn      = (esdZDC->GetZDCEMEnergy(0)+esdZDC->GetZDCEMEnergy(1))/8.;
  

  ///////////////////
  //Event Selection//
  ///////////////////
  Double_t deltaTdc = 0;
  Double_t sumTdc   = 0;
   
  for(Int_t i = 0; i < 4; ++i) {
    if(esdZDC->GetZDCTDCData(10,i) != 0) {
      Double_t  tdcCnoCorr = 0.025*(esdZDC->GetZDCTDCData(10,i)-
				    esdZDC->GetZDCTDCData(14,i)); 
      Double_t  tdcC       = esdZDC->GetZDCTDCCorrected(10,i); 
      for(Int_t j = 0; j < 4; ++j) {
	if(esdZDC->GetZDCTDCData(12,j) != 0) {
	  Double_t   tdcAnoCorr = 0.025*(esdZDC->GetZDCTDCData(12,j)-
					 esdZDC->GetZDCTDCData(14,j));
	  Double_t   tdcA       = esdZDC->GetZDCTDCCorrected(12,j);
	  if(esdZDC->TestBit(AliESDZDC::kCorrectedTDCFilled)) {	    
	    deltaTdc = tdcC-tdcA;
	    sumTdc   = tdcC+tdcA;
	  }
	  else {
	    deltaTdc = tdcCnoCorr-tdcAnoCorr;
	    sumTdc   = tdcCnoCorr+tdcAnoCorr;
	  }
	}
      }
    }
  }

  ///////////////////////
  //Global Event Filter//
  ///////////////////////
  Bool_t zdcAccSat[21];

  // Bool_t zdcAccSatRunClass[21];
  for(Int_t k = -10; k <= 10; k++) zdcAccSat[k+10] = kFALSE;


  if(deltaTdc!=0. || sumTdc!=0.) {
    for (Int_t k = -10; k <= 10; ++k) {
      Float_t zsat = 2.55F * k;
		
      if(k==0) {
	if(((deltaTdc-kZDCrefDelta)*(deltaTdc-kZDCrefDelta)/
	    (kZDCsigmaDelta*kZDCsigmaDelta)+
	    (sumTdc  -kZDCrefSum  ) * (sumTdc  -kZDCrefSum  ) / 
	    (kZDCsigmaSum  *kZDCsigmaSum))<=1.0) {
	  zdcAccSat[k+10] = kTRUE;
	}	
      }
      else {
	if(((deltaTdc-kZDCrefDelta-zsat)*(deltaTdc-kZDCrefDelta-zsat)/
	    (kZDCsigmaDeltaSat*kZDCsigmaDeltaSat)+
	    (sumTdc  -kZDCrefSum  -zsat)*(sumTdc  -kZDCrefSum  -zsat)/
	    (kZDCsigmaSumSat  *kZDCsigmaSumSat  ))<=1.0) {
	  zdcAccSat[k+10] = kTRUE;
	}		
      }
    }
  }
   
  for(Int_t k=-10;k<=10;++k) {
    if(zdcAccSat[k+10] && k!=0 ) {
      //std::cout<<"Displaced vertex at "<<37.5*(Float_t)k<<" cm"<<std::endl; 
      zvtx = 37.5*(Float_t)k;
    }
  }
 
  return zvtx;
 
}
//____________________________________________________________________
Double_t
AliDisplacedVertexSelection::CalculateDisplacedVertexCent(const AliESDEvent* esd) const 
{ 
  Float_t kZDCrefSum        =-66.5;
  Float_t kZDCrefDelta      =-2.10;
  Float_t kZDCsigmaSum      = 3.25;
  Float_t kZDCsigmaDelta    = 2.25;
  Float_t kZDCsigmaSumSat   = 2.50;
  Float_t kZDCsigmaDeltaSat = 1.35;  
  
  AliESDZDC* esdZDC = esd -> GetESDZDC();
  
  Int_t runnumber = esd->GetRunNumber();
  Double_t fcurrentL3 = esd->GetCurrentL3();
  Double_t fcurrentDipo = esd->GetCurrentDip();
  Double_t cent = -1;
  
  Double_t fzdcEn      = (esdZDC->GetZDCN1Energy()+
			  esdZDC->GetZDCP1Energy()+
			  esdZDC->GetZDCN2Energy()+
			  esdZDC->GetZDCP2Energy());
  Double_t fzemEn      = (esdZDC->GetZDCEMEnergy(0)+
			  esdZDC->GetZDCEMEnergy(1))/8.;
  

  ///////////////////
  //Event Selection//
  ///////////////////
  Double_t deltaTdc = 0;
  Double_t sumTdc   = 0;
   
  for(Int_t i = 0; i < 4; ++i) {
    if(esdZDC->GetZDCTDCData(10,i) != 0) {
      Double_t  tdcCnoCorr = 0.025*(esdZDC->GetZDCTDCData(10,i)-
				    esdZDC->GetZDCTDCData(14,i)); 
      Double_t  tdcC       = esdZDC->GetZDCTDCCorrected(10,i); 
      for(Int_t j = 0; j < 4; ++j) {
	if(esdZDC->GetZDCTDCData(12,j) != 0) {
	  Double_t   tdcAnoCorr = 0.025*(esdZDC->GetZDCTDCData(12,j)-
					 esdZDC->GetZDCTDCData(14,j));
	  Double_t   tdcA       = esdZDC->GetZDCTDCCorrected(12,j);
	  if(esdZDC->TestBit(AliESDZDC::kCorrectedTDCFilled)) {	    
	    deltaTdc = tdcC-tdcA;
	    sumTdc   = tdcC+tdcA;
	  }
	  else {
	    deltaTdc = tdcCnoCorr-tdcAnoCorr;
	    sumTdc   = tdcCnoCorr+tdcAnoCorr;
	  }
	}
      }
    }
  }

  ///////////////////////
  //Global Event Filter//
  ///////////////////////
  Bool_t zdcAccSat[21];
  // Bool_t zdcAccSatRunClass[21];
  for(Int_t k = -10; k <= 10; k++)  zdcAccSat[k+10]         = kFALSE;

  if(deltaTdc!=0. || sumTdc!=0.) {
    for(Int_t k = -10; k <= 10; ++k) {
      Float_t zsat = 2.55*(Float_t)k;
		
      if(k==0) {
	if(((deltaTdc-kZDCrefDelta)*(deltaTdc-kZDCrefDelta)/
	    (kZDCsigmaDelta*kZDCsigmaDelta)+
	    (sumTdc  -kZDCrefSum  )*(sumTdc  -kZDCrefSum  )/
	    (kZDCsigmaSum  *kZDCsigmaSum))<=1.0) {
	  zdcAccSat[k+10] = kTRUE;
	}	
      }
      else {
	if(((deltaTdc-kZDCrefDelta-zsat)*(deltaTdc-kZDCrefDelta-zsat)/
	    (kZDCsigmaDeltaSat*kZDCsigmaDeltaSat)+
	    (sumTdc  -kZDCrefSum  -zsat)*(sumTdc  -kZDCrefSum  -zsat)/
	    (kZDCsigmaSumSat  *kZDCsigmaSumSat  ))<=1.0) {
	  zdcAccSat[k+10] = kTRUE;
	}		
	
      }
      
    }
  }

  Float_t kZEMcorrPlusPlus[21]   = {0.6840,0.7879,0.8722,0.9370,0.9837,1.0137,
				    1.0292,1.0327,1.0271,1.0152,1.0000,0.9844,
				    0.9714,0.9634,0.9626,0.9708,0.9891,1.0181,
				    1.0574,1.1060,1.1617};
  Float_t kZEMcorrMoinsMoins[21] = {0.7082,0.8012,0.8809,0.9447,0.9916,1.0220,
				    1.0372,1.0395,1.0318,1.0174,1.0000,0.9830,
				    0.9699,0.9635,0.9662,0.9794,1.0034,1.0371,
				    1.0782,1.1224,1.1634};
  ///////////////////
  //ZemEn correction//
  ////////////////////

  for(Int_t k = -10; k <= 10; ++k) {        
    if(zdcAccSat[k+10]) { 
      //std::cout<<"Vertex at "<<37.5*(Float_t)k<<"cm from cent"<<std::endl;
      if(fcurrentDipo>0 && fcurrentL3>0) {
	fzemEn /= kZEMcorrPlusPlus[k+10];
      }
      if(fcurrentDipo<0 && fcurrentL3<0) {
	fzemEn /= kZEMcorrMoinsMoins[k+10];
      }
    }          
  }
 
  ////////////////////////
  //Centrality selection//
  ////////////////////////
  Double_t fzdcPercentile = -1;
  Float_t slope;
 
 
  if(runnumber < 137722 && runnumber >= 137161 ) {
    if(fzemEn > 345.) {
      slope = (fzdcEn + 16992.)/(fzemEn - 345.);
      slope += 2.23902e+02;
      fzdcPercentile = (TMath::ATan(slope) - 1.56667)/9.49434e-05;
      if (fzdcPercentile<0) fzdcPercentile = 0;
    }
    else fzdcPercentile = 100;                }
  else if(runnumber < 139172 && runnumber >= 137722) {
    if(fzemEn > 295.) {
      slope = (fzdcEn + 15000.)/(fzemEn - 295.);
      slope += 2.23660e+02;
      fzdcPercentile = (TMath::ATan(slope) - 1.56664)/8.99571e-05;
      if (fzdcPercentile<0) fzdcPercentile = 0;
    }
    else fzdcPercentile = 100;          
  }
  else if(runnumber >= 139172) {
    if(fzemEn > 345.) {
      slope = (fzdcEn + 16992.)/(fzemEn - 345.);
      slope += 2.04789e+02;
      fzdcPercentile = (TMath::ATan(slope) - 1.56629)/1.02768e-04;
      if (fzdcPercentile<0) fzdcPercentile = 0;
    }                               
    else fzdcPercentile = 100;           
  }
  else {
    if(fzemEn > 345.) {
      slope = (fzdcEn + 16992.)/(fzemEn - 345.);
      slope += 2.04789e+02;
      fzdcPercentile = (TMath::ATan(slope) - 1.56629)/1.02768e-04;
      if(fzdcPercentile<0) fzdcPercentile = 0;
    }
    else fzdcPercentile = 100;                  
  }
  
  if(fzdcPercentile > 0 && fzdcPercentile <100) {
    //std::cout<<fzdcPercentile<<std::endl;
    cent = fzdcPercentile;
  }
  return cent;
}
//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
65abd48b	1	#include "AliDisplacedVertexSelection.h"
	2	#include <iostream>
	3	#include <TROOT.h>
	4	#include "AliESDEvent.h"
	5	#include "AliESDZDC.h"
	6	ClassImp(AliDisplacedVertexSelection)
	7	#if 0
	8	; // This is for Emacs
	9	#endif
	10
	11	//____________________________________________________________________
	12	AliDisplacedVertexSelection::AliDisplacedVertexSelection()
241cca4d	13	: TObject(),
	14	fVertexZ(9999),
	15	fCent(100)
65abd48b	16	{
	17	}
	18	//____________________________________________________________________
	19	AliDisplacedVertexSelection::AliDisplacedVertexSelection(const AliDisplacedVertexSelection& o)
241cca4d	20	: TObject(o),
	21	fVertexZ(9999),
	22	fCent(100)
65abd48b	23	{
	24	}
	25	//____________________________________________________________________
	26	AliDisplacedVertexSelection&
	27	AliDisplacedVertexSelection::operator=(const AliDisplacedVertexSelection& o)
	28	{
	29	if (&o == this) return *this;
	30	return *this;
	31	}
	32
	33	//____________________________________________________________________
	34	void
	35	AliDisplacedVertexSelection::Output(TList* /l/, const char* /* name*/) const
	36	{
	37	}
	38
	39	//____________________________________________________________________
	40	void
	41	AliDisplacedVertexSelection::Print(Option_t*) const
	42	{
6ab100ec	43	#if 0
65abd48b	44	char ind[gROOT->GetDirLevel()+1];
	45	for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
	46	ind[gROOT->GetDirLevel()] = '\0';
	47	std::cout << std::boolalpha
	48	<< std::noboolalpha << std::endl;
6ab100ec	49	#endif
65abd48b	50	}
241cca4d	51
65abd48b	52	//____________________________________________________________________
241cca4d	53	Bool_t
	54	AliDisplacedVertexSelection::Process(const AliESDEvent* esd)
	55	{
	56	fVertexZ = 9999; // Default vertex value
	57	fCent = 100; // Default centrality value
	58
	59	// Some constants
	60	const Float_t kZDCrefSum = -66.5;
	61	const Float_t kZDCrefDelta = -2.10;
	62	const Float_t kZDCsigmaSum = 3.25;
	63	const Float_t kZDCsigmaDelta = 2.25;
	64	const Float_t kZDCsigmaSumSat = 2.50;
	65	const Float_t kZDCsigmaDeltaSat = 1.35;
	66
	67	// Corrections for magnetic fields
	68	const Float_t kZEMcorrPlusPlus[21] = {0.6840,0.7879,0.8722,
	69	0.9370,0.9837,1.0137,
	70	1.0292,1.0327,1.0271,
	71	1.0152,1.0000,0.9844,
	72	0.9714,0.9634,0.9626,
	73	0.9708,0.9891,1.0181,
	74	1.0574,1.1060,1.1617};
	75	const Float_t kZEMcorrMoinsMoins[21] = {0.7082,0.8012,0.8809,
	76	0.9447,0.9916,1.0220,
	77	1.0372,1.0395,1.0318,
	78	1.0174,1.0000,0.9830,
	79	0.9699,0.9635,0.9662,
	80	0.9794,1.0034,1.0371,
	81	1.0782,1.1224,1.1634};
	82
	83	// --- Get the ESD object ------------------------------------------
	84	AliESDZDC* esdZDC = esd->GetESDZDC();
	85	if (!esdZDC) {
	86	AliWarning("No ZDC ESD object!");
	87	return false;
	88	}
	89	Double_t currentL3 = esd->GetCurrentL3();
	90	Double_t currentDipo = esd->GetCurrentDip();
	91
	92	// --- ZDC and ZEM energy signal -----------------------------------
	93	Double_t zdcEn = (esdZDC->GetZDCN1Energy()+
	94	esdZDC->GetZDCP1Energy()+
	95	esdZDC->GetZDCN2Energy()+
	96	esdZDC->GetZDCP2Energy());
	97	Double_t zemEn = (esdZDC->GetZDCEMEnergy(0)+
	98	esdZDC->GetZDCEMEnergy(1))/8.;
	99
	100	// --- Calculate DeltaT and sumT -----------------------------------
	101	Double_t deltaTdc = 0;
	102	Double_t sumTdc = 0;
	103
	104	for(Int_t i = 0; i < 4; ++i) {
	105	if(esdZDC->GetZDCTDCData(10,i) != 0) {
	106	Double_t tdcCnoCorr = 0.025*(esdZDC->GetZDCTDCData(10,i)-
	107	esdZDC->GetZDCTDCData(14,i));
	108	Double_t tdcC = esdZDC->GetZDCTDCCorrected(10,i);
	109	for(Int_t j = 0; j < 4; ++j) {
	110	if(esdZDC->GetZDCTDCData(12,j) != 0) {
	111	Double_t tdcAnoCorr = 0.025*(esdZDC->GetZDCTDCData(12,j)-
	112	esdZDC->GetZDCTDCData(14,j));
	113	Double_t tdcA = esdZDC->GetZDCTDCCorrected(12,j);
	114	if(esdZDC->TestBit(AliESDZDC::kCorrectedTDCFilled)) {
	115	deltaTdc = tdcC-tdcA;
	116	sumTdc = tdcC+tdcA;
117	}
118	else {
119	deltaTdc = tdcCnoCorr-tdcAnoCorr;
120	sumTdc = tdcCnoCorr+tdcAnoCorr;
121	}
122	}
123	}
124	}
125	}
126
127	// --- Find the vertex ---------------------------------------------
128	if(deltaTdc!=0. \|\| sumTdc!=0.) {
129	for (Int_t k = -10; k <= 10; ++k) {
130	Float_t zsat = 2.55F * k;
131	Float_t delta = (k == 0 ? kZDCsigmaDelta : kZDCsigmaDeltaSat);
132	Float_t sum = (k == 0 ? kZDCsigmaSum : kZDCsigmaSumSat);
133	Float_t dT = deltaTdc - kZDCrefDelta - zsat;
134	Float_t sT = sumTdc - kZDCrefSum - zsat;
135	Float_t check = dT * dT / delta / delta + sT * sT / sum / sum;
136	if (check > 1.0 \|\| k == 0) continue;
137
138	// Set the vertex
139	fVertexZ = 37.5 * k;
140
141	// Correct zem energy
142	if(currentDipo>0 && currentL3>0) zemEn /= kZEMcorrPlusPlus[k+10];
143	if(currentDipo<0 && currentL3<0) zemEn /= kZEMcorrMoinsMoins[k+10];
144	}
145	}
146
147	// --- Calculate the centrality ------------------------------------
148	Float_t c1, c2, c3, c4, c5;
149	Int_t runnumber = esd->GetRunNumber();
150	if (runnumber < 137722 && runnumber >= 137161 ) {
151	c1 = 16992;
152	c2 = 345;
153	c3 = 2.23902e+02;
154	c4 = 1.56667;
155	c5 = 9.49434e-05;
156	}
157	else if (runnumber < 139172 && runnumber >= 137722) {
158	c1 = 15000;
159	c2 = 295;
160	c3 = 2.23660e+02;
161	c4 = 1.56664;
162	c5 = 8.99571e-05;
163	}
164	else if (runnumber >= 139172) {
165	c1 = 16992;
166	c2 = 345;
167	c3 = 2.04789e+02;
168	c4 = 1.56629;
169	c5 = 1.02768e-04;
170	}
171	else {
172	c1 = 16992;
173	c2 = 345;
174	c3 = 2.04789e+02;
175	c4 = 1.56629;
176	c5 = 1.02768e-04;
177	}
178
179	if (zemEn > c2) {
180	Float_t slope = (zdcEn + c1) / (zemEn - c2) + c3;
181	Float_t zdcCent = (TMath::ATan(slope) - c4) / c5;
182	if (zdcCent >= 0) fCent = zdcCent;
183	}
65abd48b	184
241cca4d	185	return true;
	186	}
	187
	188	//____________________________________________________________________
65abd48b	189	Double_t
	190	AliDisplacedVertexSelection::CheckDisplacedVertex(const AliESDEvent* esd) const {
	191	// Code taken directly from M.Guilbaud.
241cca4d	192	AliWarning("This method is deprecated");
65abd48b	193	Double_t zvtx = 9999.;
	194
	195	Float_t kZDCrefSum =-66.5;
	196	Float_t kZDCrefDelta =-2.10;
	197	Float_t kZDCsigmaSum = 3.25;
	198	Float_t kZDCsigmaDelta = 2.25;
	199	Float_t kZDCsigmaSumSat = 2.50;
	200	Float_t kZDCsigmaDeltaSat = 1.35;
	201
	202	AliESDZDC* esdZDC = esd -> GetESDZDC();
	203
	204	/* Double_t zdcNCEnergy = esdZDC->GetZDCN1Energy();
	205	Double_t fZpcEnergy = esdZDC->GetZDCP1Energy();
	206	Double_t fZnaEnergy = esdZDC->GetZDCN2Energy();
	207	Double_t fZpaEnergy = esdZDC->GetZDCP2Energy();
	208	Double_t fZem1Energy = esdZDC->GetZDCEMEnergy(0);
	209	Double_t fZem2Energy = esdZDC->GetZDCEMEnergy(1);*/
	210
	211	//Double_t fzdcEn = (esdZDC->GetZDCN1Energy()+esdZDC->GetZDCP1Energy()+esdZDC->GetZDCN2Energy()+esdZDC->GetZDCP2Energy());
	212	//Double_t fzemEn = (esdZDC->GetZDCEMEnergy(0)+esdZDC->GetZDCEMEnergy(1))/8.;
	213
	214
241cca4d	215
	216	///////////////////
	217	//Event Selection//
	218	///////////////////
	219	Double_t deltaTdc = 0;
	220	Double_t sumTdc = 0;
65abd48b	221
241cca4d	222	for(Int_t i = 0; i < 4; ++i) {
	223	if(esdZDC->GetZDCTDCData(10,i) != 0) {
	224	Double_t tdcCnoCorr = 0.025*(esdZDC->GetZDCTDCData(10,i)-
	225	esdZDC->GetZDCTDCData(14,i));
	226	Double_t tdcC = esdZDC->GetZDCTDCCorrected(10,i);
	227	for(Int_t j = 0; j < 4; ++j) {
	228	if(esdZDC->GetZDCTDCData(12,j) != 0) {
	229	Double_t tdcAnoCorr = 0.025*(esdZDC->GetZDCTDCData(12,j)-
	230	esdZDC->GetZDCTDCData(14,j));
	231	Double_t tdcA = esdZDC->GetZDCTDCCorrected(12,j);
	232	if(esdZDC->TestBit(AliESDZDC::kCorrectedTDCFilled)) {
	233	deltaTdc = tdcC-tdcA;
	234	sumTdc = tdcC+tdcA;
	235	}
	236	else {
	237	deltaTdc = tdcCnoCorr-tdcAnoCorr;
	238	sumTdc = tdcCnoCorr+tdcAnoCorr;
	239	}
	240	}
	241	}
	242	}
	243	}
65abd48b	244
241cca4d	245	///////////////////////
	246	//Global Event Filter//
	247	///////////////////////
	248	Bool_t zdcAccSat[21];
65abd48b	249
241cca4d	250	// Bool_t zdcAccSatRunClass[21];
241cca4d	251	for(Int_t k = -10; k <= 10; k++) zdcAccSat[k+10] = kFALSE;
65abd48b	252
241cca4d	253
	254	if(deltaTdc!=0. \|\| sumTdc!=0.) {
	255	for (Int_t k = -10; k <= 10; ++k) {
	256	Float_t zsat = 2.55F * k;
65abd48b	257
241cca4d	258	if(k==0) {
	259	if(((deltaTdc-kZDCrefDelta)*(deltaTdc-kZDCrefDelta)/
	260	(kZDCsigmaDelta*kZDCsigmaDelta)+
	261	(sumTdc -kZDCrefSum ) * (sumTdc -kZDCrefSum ) /
	262	(kZDCsigmaSum *kZDCsigmaSum))<=1.0) {
	263	zdcAccSat[k+10] = kTRUE;
	264	}
	265	}
	266	else {
	267	if(((deltaTdc-kZDCrefDelta-zsat)*(deltaTdc-kZDCrefDelta-zsat)/
	268	(kZDCsigmaDeltaSat*kZDCsigmaDeltaSat)+
	269	(sumTdc -kZDCrefSum -zsat)*(sumTdc -kZDCrefSum -zsat)/
	270	(kZDCsigmaSumSat *kZDCsigmaSumSat ))<=1.0) {
	271	zdcAccSat[k+10] = kTRUE;
	272	}
	273	}
	274	}
65abd48b	275	}
65abd48b	276
241cca4d	277	for(Int_t k=-10;k<=10;++k) {
	278	if(zdcAccSat[k+10] && k!=0 ) {
	279	//std::cout<<"Displaced vertex at "<<37.5*(Float_t)k<<" cm"<<std::endl;
	280	zvtx = 37.5*(Float_t)k;
	281	}
	282	}
65abd48b	283
241cca4d	284	return zvtx;
65abd48b	285
	286	}
	287	//____________________________________________________________________
	288	Double_t
	289	AliDisplacedVertexSelection::CalculateDisplacedVertexCent(const AliESDEvent* esd) const
	290	{
	291	Float_t kZDCrefSum =-66.5;
	292	Float_t kZDCrefDelta =-2.10;
	293	Float_t kZDCsigmaSum = 3.25;
	294	Float_t kZDCsigmaDelta = 2.25;
	295	Float_t kZDCsigmaSumSat = 2.50;
	296	Float_t kZDCsigmaDeltaSat = 1.35;
	297
	298	AliESDZDC* esdZDC = esd -> GetESDZDC();
	299
	300	Int_t runnumber = esd->GetRunNumber();
	301	Double_t fcurrentL3 = esd->GetCurrentL3();
	302	Double_t fcurrentDipo = esd->GetCurrentDip();
	303	Double_t cent = -1;
	304
241cca4d	305	Double_t fzdcEn = (esdZDC->GetZDCN1Energy()+
	306	esdZDC->GetZDCP1Energy()+
	307	esdZDC->GetZDCN2Energy()+
	308	esdZDC->GetZDCP2Energy());
	309	Double_t fzemEn = (esdZDC->GetZDCEMEnergy(0)+
	310	esdZDC->GetZDCEMEnergy(1))/8.;
65abd48b	311
65abd48b	312
241cca4d	313	///////////////////
	314	//Event Selection//
	315	///////////////////
	316	Double_t deltaTdc = 0;
	317	Double_t sumTdc = 0;
65abd48b	318
241cca4d	319	for(Int_t i = 0; i < 4; ++i) {
	320	if(esdZDC->GetZDCTDCData(10,i) != 0) {
	321	Double_t tdcCnoCorr = 0.025*(esdZDC->GetZDCTDCData(10,i)-
	322	esdZDC->GetZDCTDCData(14,i));
	323	Double_t tdcC = esdZDC->GetZDCTDCCorrected(10,i);
	324	for(Int_t j = 0; j < 4; ++j) {
	325	if(esdZDC->GetZDCTDCData(12,j) != 0) {
	326	Double_t tdcAnoCorr = 0.025*(esdZDC->GetZDCTDCData(12,j)-
	327	esdZDC->GetZDCTDCData(14,j));
	328	Double_t tdcA = esdZDC->GetZDCTDCCorrected(12,j);
	329	if(esdZDC->TestBit(AliESDZDC::kCorrectedTDCFilled)) {
	330	deltaTdc = tdcC-tdcA;
	331	sumTdc = tdcC+tdcA;
	332	}
	333	else {
	334	deltaTdc = tdcCnoCorr-tdcAnoCorr;
	335	sumTdc = tdcCnoCorr+tdcAnoCorr;
	336	}
	337	}
	338	}
	339	}
	340	}
65abd48b	341
241cca4d	342	///////////////////////
	343	//Global Event Filter//
	344	///////////////////////
	345	Bool_t zdcAccSat[21];
	346	// Bool_t zdcAccSatRunClass[21];
	347	for(Int_t k = -10; k <= 10; k++) zdcAccSat[k+10] = kFALSE;
65abd48b	348
241cca4d	349	if(deltaTdc!=0. \|\| sumTdc!=0.) {
	350	for(Int_t k = -10; k <= 10; ++k) {
	351	Float_t zsat = 2.55*(Float_t)k;
65abd48b	352
241cca4d	353	if(k==0) {
	354	if(((deltaTdc-kZDCrefDelta)*(deltaTdc-kZDCrefDelta)/
	355	(kZDCsigmaDelta*kZDCsigmaDelta)+
	356	(sumTdc -kZDCrefSum )*(sumTdc -kZDCrefSum )/
	357	(kZDCsigmaSum *kZDCsigmaSum))<=1.0) {
	358	zdcAccSat[k+10] = kTRUE;
	359	}
	360	}
	361	else {
	362	if(((deltaTdc-kZDCrefDelta-zsat)*(deltaTdc-kZDCrefDelta-zsat)/
	363	(kZDCsigmaDeltaSat*kZDCsigmaDeltaSat)+
	364	(sumTdc -kZDCrefSum -zsat)*(sumTdc -kZDCrefSum -zsat)/
	365	(kZDCsigmaSumSat *kZDCsigmaSumSat ))<=1.0) {
	366	zdcAccSat[k+10] = kTRUE;
	367	}
	368
	369	}
	370
	371	}
65abd48b	372	}
65abd48b	373
241cca4d	374	Float_t kZEMcorrPlusPlus[21] = {0.6840,0.7879,0.8722,0.9370,0.9837,1.0137,
	375	1.0292,1.0327,1.0271,1.0152,1.0000,0.9844,
	376	0.9714,0.9634,0.9626,0.9708,0.9891,1.0181,
	377	1.0574,1.1060,1.1617};
	378	Float_t kZEMcorrMoinsMoins[21] = {0.7082,0.8012,0.8809,0.9447,0.9916,1.0220,
	379	1.0372,1.0395,1.0318,1.0174,1.0000,0.9830,
	380	0.9699,0.9635,0.9662,0.9794,1.0034,1.0371,
	381	1.0782,1.1224,1.1634};
	382	///////////////////
	383	//ZemEn correction//
	384	////////////////////
65abd48b	385
241cca4d	386	for(Int_t k = -10; k <= 10; ++k) {
	387	if(zdcAccSat[k+10]) {
	388	//std::cout<<"Vertex at "<<37.5*(Float_t)k<<"cm from cent"<<std::endl;
	389	if(fcurrentDipo>0 && fcurrentL3>0) {
	390	fzemEn /= kZEMcorrPlusPlus[k+10];
	391	}
	392	if(fcurrentDipo<0 && fcurrentL3<0) {
	393	fzemEn /= kZEMcorrMoinsMoins[k+10];
	394	}
	395	}
	396	}
65abd48b	397
241cca4d	398	////////////////////////
	399	//Centrality selection//
	400	////////////////////////
	401	Double_t fzdcPercentile = -1;
	402	Float_t slope;
65abd48b	403
65abd48b	404
241cca4d	405	if(runnumber < 137722 && runnumber >= 137161 ) {
	406	if(fzemEn > 345.) {
	407	slope = (fzdcEn + 16992.)/(fzemEn - 345.);
	408	slope += 2.23902e+02;
	409	fzdcPercentile = (TMath::ATan(slope) - 1.56667)/9.49434e-05;
	410	if (fzdcPercentile<0) fzdcPercentile = 0;
	411	}
	412	else fzdcPercentile = 100; }
	413	else if(runnumber < 139172 && runnumber >= 137722) {
	414	if(fzemEn > 295.) {
	415	slope = (fzdcEn + 15000.)/(fzemEn - 295.);
	416	slope += 2.23660e+02;
	417	fzdcPercentile = (TMath::ATan(slope) - 1.56664)/8.99571e-05;
	418	if (fzdcPercentile<0) fzdcPercentile = 0;
	419	}
	420	else fzdcPercentile = 100;
	421	}
	422	else if(runnumber >= 139172) {
	423	if(fzemEn > 345.) {
	424	slope = (fzdcEn + 16992.)/(fzemEn - 345.);
	425	slope += 2.04789e+02;
	426	fzdcPercentile = (TMath::ATan(slope) - 1.56629)/1.02768e-04;
	427	if (fzdcPercentile<0) fzdcPercentile = 0;
	428	}
	429	else fzdcPercentile = 100;
	430	}
	431	else {
	432	if(fzemEn > 345.) {
	433	slope = (fzdcEn + 16992.)/(fzemEn - 345.);
	434	slope += 2.04789e+02;
	435	fzdcPercentile = (TMath::ATan(slope) - 1.56629)/1.02768e-04;
	436	if(fzdcPercentile<0) fzdcPercentile = 0;
	437	}
	438	else fzdcPercentile = 100;
	439	}
	440
	441	if(fzdcPercentile > 0 && fzdcPercentile <100) {
	442	//std::cout<<fzdcPercentile<<std::endl;
	443	cent = fzdcPercentile;
	444	}
	445	return cent;
65abd48b	446	}
	447	//____________________________________________________________________
	448	//
	449	// EOF
	450	//