[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()
{
}
//____________________________________________________________________
AliDisplacedVertexSelection::AliDisplacedVertexSelection(const AliDisplacedVertexSelection& o)
  : TObject(o)
{
}
//____________________________________________________________________
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
}
//____________________________________________________________________

Double_t
AliDisplacedVertexSelection::CheckDisplacedVertex(const AliESDEvent* esd) const {
  // Code taken directly from M.Guilbaud.
  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;
				}
	                 }
	          }
          }
   }
//if(TMath::Abs(deltaTdc) > 0)
//std::cout<<deltaTdc<<"   "<<sumTdc<<std::endl;
///////////////////////
//Global Event Filter//
///////////////////////
 Bool_t zdcAccSat[21];
 // Bool_t zdcAccSatRunClass[21];
 for(Int_t k = -10; k <= 10; k++) 
   {
      zdcAccSat[k+10]         = kFALSE;
      // zdcAccSatRunClass[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;
		       }		
	  
	       }
	  
         }
  }
   
 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 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;
				}
	                 }
	          }
          }
   }
//if(TMath::Abs(deltaTdc) > 0)
//std::cout<<deltaTdc<<"   "<<sumTdc<<std::endl;
///////////////////////
//Global Event Filter//
///////////////////////
 Bool_t zdcAccSat[21];
 // Bool_t zdcAccSatRunClass[21];
for(Int_t k = -10; k <= 10; k++) 
   {
      zdcAccSat[k+10]         = kFALSE;
      // zdcAccSatRunClass[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()
	13	: TObject()
	14	{
	15	}
	16	//____________________________________________________________________
	17	AliDisplacedVertexSelection::AliDisplacedVertexSelection(const AliDisplacedVertexSelection& o)
	18	: TObject(o)
	19	{
	20	}
	21	//____________________________________________________________________
	22	AliDisplacedVertexSelection&
	23	AliDisplacedVertexSelection::operator=(const AliDisplacedVertexSelection& o)
	24	{
	25	if (&o == this) return *this;
	26	return *this;
	27	}
	28
	29	//____________________________________________________________________
	30	void
	31	AliDisplacedVertexSelection::Output(TList* /l/, const char* /* name*/) const
	32	{
	33	}
	34
	35	//____________________________________________________________________
	36	void
	37	AliDisplacedVertexSelection::Print(Option_t*) const
	38	{
6ab100ec	39	#if 0
65abd48b	40	char ind[gROOT->GetDirLevel()+1];
	41	for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
	42	ind[gROOT->GetDirLevel()] = '\0';
	43	std::cout << std::boolalpha
	44	<< std::noboolalpha << std::endl;
6ab100ec	45	#endif
65abd48b	46	}
	47	//____________________________________________________________________
	48
	49	Double_t
	50	AliDisplacedVertexSelection::CheckDisplacedVertex(const AliESDEvent* esd) const {
	51	// Code taken directly from M.Guilbaud.
	52	Double_t zvtx = 9999.;
	53
	54	Float_t kZDCrefSum =-66.5;
	55	Float_t kZDCrefDelta =-2.10;
	56	Float_t kZDCsigmaSum = 3.25;
	57	Float_t kZDCsigmaDelta = 2.25;
	58	Float_t kZDCsigmaSumSat = 2.50;
	59	Float_t kZDCsigmaDeltaSat = 1.35;
	60
	61	AliESDZDC* esdZDC = esd -> GetESDZDC();
	62
	63	/* Double_t zdcNCEnergy = esdZDC->GetZDCN1Energy();
	64	Double_t fZpcEnergy = esdZDC->GetZDCP1Energy();
	65	Double_t fZnaEnergy = esdZDC->GetZDCN2Energy();
	66	Double_t fZpaEnergy = esdZDC->GetZDCP2Energy();
	67	Double_t fZem1Energy = esdZDC->GetZDCEMEnergy(0);
	68	Double_t fZem2Energy = esdZDC->GetZDCEMEnergy(1);*/
	69
	70	//Double_t fzdcEn = (esdZDC->GetZDCN1Energy()+esdZDC->GetZDCP1Energy()+esdZDC->GetZDCN2Energy()+esdZDC->GetZDCP2Energy());
	71	//Double_t fzemEn = (esdZDC->GetZDCEMEnergy(0)+esdZDC->GetZDCEMEnergy(1))/8.;
	72
	73
	74	///////////////////
	75	//Event Selection//
	76	///////////////////
	77	Double_t deltaTdc = 0;
	78	Double_t sumTdc = 0;
	79
	80	for(Int_t i = 0; i < 4; ++i)
	81	{
	82	if(esdZDC->GetZDCTDCData(10,i) != 0)
	83	{
	84	Double_t tdcCnoCorr = 0.025*(esdZDC->GetZDCTDCData(10,i)-esdZDC->GetZDCTDCData(14,i));
	85	Double_t tdcC = esdZDC->GetZDCTDCCorrected(10,i);
	86	for(Int_t j = 0; j < 4; ++j)
	87	{
	88	if(esdZDC->GetZDCTDCData(12,j) != 0)
	89	{
	90	Double_t tdcAnoCorr = 0.025*(esdZDC->GetZDCTDCData(12,j)-esdZDC->GetZDCTDCData(14,j));
	91	Double_t tdcA = esdZDC->GetZDCTDCCorrected(12,j);
	92	if(esdZDC->TestBit(AliESDZDC::kCorrectedTDCFilled))
	93	{
	94	deltaTdc = tdcC-tdcA;
	95	sumTdc = tdcC+tdcA;
	96	}
	97	else
	98	{
	99	deltaTdc = tdcCnoCorr-tdcAnoCorr;
	100	sumTdc = tdcCnoCorr+tdcAnoCorr;
	101	}
	102	}
	103	}
	104	}
	105	}
	106	//if(TMath::Abs(deltaTdc) > 0)
	107	//std::cout<<deltaTdc<<" "<<sumTdc<<std::endl;
	108	///////////////////////
	109	//Global Event Filter//
110	///////////////////////
111	Bool_t zdcAccSat[21];
6ab100ec	112	// Bool_t zdcAccSatRunClass[21];
6ab100ec	113	for(Int_t k = -10; k <= 10; k++)
65abd48b	114	{
65abd48b	115	zdcAccSat[k+10] = kFALSE;
6ab100ec	116	// zdcAccSatRunClass[k+10] = kFALSE;
65abd48b	117	}
	118
	119
	120	if(deltaTdc!=0. \|\| sumTdc!=0.)
	121	{
	122
	123	for(Int_t k = -10; k <= 10; ++k)
	124	{
	125	Float_t zsat = 2.55*(Float_t)k;
	126
	127	if(k==0)
	128	{
	129	if(((deltaTdc-kZDCrefDelta)(deltaTdc-kZDCrefDelta)/(kZDCsigmaDeltakZDCsigmaDelta)+
	130	(sumTdc -kZDCrefSum )(sumTdc -kZDCrefSum )/(kZDCsigmaSum kZDCsigmaSum))<=1.0)
	131	{
	132	zdcAccSat[k+10] = kTRUE;
	133	}
	134	}
	135	else
	136	{
	137	if(((deltaTdc-kZDCrefDelta-zsat)(deltaTdc-kZDCrefDelta-zsat)/(kZDCsigmaDeltaSatkZDCsigmaDeltaSat)+
	138	(sumTdc -kZDCrefSum -zsat)(sumTdc -kZDCrefSum -zsat)/(kZDCsigmaSumSat kZDCsigmaSumSat ))<=1.0)
	139	{
	140	zdcAccSat[k+10] = kTRUE;
	141	}
	142
	143	}
	144
	145	}
	146	}
	147
	148	for(Int_t k=-10;k<=10;++k) {
	149	if(zdcAccSat[k+10] && k!=0 ) {
	150	//std::cout<<"Displaced vertex at "<<37.5*(Float_t)k<<" cm"<<std::endl;
	151	zvtx = 37.5*(Float_t)k;
	152	}
	153	}
	154
	155	return zvtx;
	156
	157	}
	158	//____________________________________________________________________
	159	Double_t
	160	AliDisplacedVertexSelection::CalculateDisplacedVertexCent(const AliESDEvent* esd) const
	161	{
	162	Float_t kZDCrefSum =-66.5;
	163	Float_t kZDCrefDelta =-2.10;
	164	Float_t kZDCsigmaSum = 3.25;
	165	Float_t kZDCsigmaDelta = 2.25;
	166	Float_t kZDCsigmaSumSat = 2.50;
	167	Float_t kZDCsigmaDeltaSat = 1.35;
	168
	169	AliESDZDC* esdZDC = esd -> GetESDZDC();
	170
	171	Int_t runnumber = esd->GetRunNumber();
	172	Double_t fcurrentL3 = esd->GetCurrentL3();
	173	Double_t fcurrentDipo = esd->GetCurrentDip();
	174	Double_t cent = -1;
	175
	176	/*Double_t zdcNCEnergy = esdZDC->GetZDCN1Energy();
	177	Double_t fZpcEnergy = esdZDC->GetZDCP1Energy();
	178	Double_t fZnaEnergy = esdZDC->GetZDCN2Energy();
	179	Double_t fZpaEnergy = esdZDC->GetZDCP2Energy();
	180	Double_t fZem1Energy = esdZDC->GetZDCEMEnergy(0);
181	Double_t fZem2Energy = esdZDC->GetZDCEMEnergy(1);*/
182
183	Double_t fzdcEn = (esdZDC->GetZDCN1Energy()+esdZDC->GetZDCP1Energy()+esdZDC->GetZDCN2Energy()+esdZDC->GetZDCP2Energy());
184	Double_t fzemEn = (esdZDC->GetZDCEMEnergy(0)+esdZDC->GetZDCEMEnergy(1))/8.;
185
186
187	///////////////////
188	//Event Selection//
189	///////////////////
190	Double_t deltaTdc = 0;
191	Double_t sumTdc = 0;
192
193	for(Int_t i = 0; i < 4; ++i)
194	{
195	if(esdZDC->GetZDCTDCData(10,i) != 0)
196	{
197	Double_t tdcCnoCorr = 0.025*(esdZDC->GetZDCTDCData(10,i)-esdZDC->GetZDCTDCData(14,i));
198	Double_t tdcC = esdZDC->GetZDCTDCCorrected(10,i);
199	for(Int_t j = 0; j < 4; ++j)
200	{
201	if(esdZDC->GetZDCTDCData(12,j) != 0)
202	{
203	Double_t tdcAnoCorr = 0.025*(esdZDC->GetZDCTDCData(12,j)-esdZDC->GetZDCTDCData(14,j));
204	Double_t tdcA = esdZDC->GetZDCTDCCorrected(12,j);
205	if(esdZDC->TestBit(AliESDZDC::kCorrectedTDCFilled))
206	{
207	deltaTdc = tdcC-tdcA;
208	sumTdc = tdcC+tdcA;
209	}
210	else
211	{
212	deltaTdc = tdcCnoCorr-tdcAnoCorr;
213	sumTdc = tdcCnoCorr+tdcAnoCorr;
214	}
215	}
216	}
217	}
218	}
219	//if(TMath::Abs(deltaTdc) > 0)
220	//std::cout<<deltaTdc<<" "<<sumTdc<<std::endl;
221	///////////////////////
222	//Global Event Filter//
223	///////////////////////
224	Bool_t zdcAccSat[21];
6ab100ec	225	// Bool_t zdcAccSatRunClass[21];
65abd48b	226	for(Int_t k = -10; k <= 10; k++)
	227	{
	228	zdcAccSat[k+10] = kFALSE;
6ab100ec	229	// zdcAccSatRunClass[k+10] = kFALSE;
65abd48b	230	}
	231
	232
	233	if(deltaTdc!=0. \|\| sumTdc!=0.)
	234	{
	235
	236	for(Int_t k = -10; k <= 10; ++k)
	237	{
	238	Float_t zsat = 2.55*(Float_t)k;
	239
	240	if(k==0)
	241	{
	242	if(((deltaTdc-kZDCrefDelta)(deltaTdc-kZDCrefDelta)/(kZDCsigmaDeltakZDCsigmaDelta)+
	243	(sumTdc -kZDCrefSum )(sumTdc -kZDCrefSum )/(kZDCsigmaSum kZDCsigmaSum))<=1.0)
	244	{
	245	zdcAccSat[k+10] = kTRUE;
	246	}
	247	}
	248	else
	249	{
	250	if(((deltaTdc-kZDCrefDelta-zsat)(deltaTdc-kZDCrefDelta-zsat)/(kZDCsigmaDeltaSatkZDCsigmaDeltaSat)+
	251	(sumTdc -kZDCrefSum -zsat)(sumTdc -kZDCrefSum -zsat)/(kZDCsigmaSumSat kZDCsigmaSumSat ))<=1.0)
	252	{
	253	zdcAccSat[k+10] = kTRUE;
	254	}
	255
	256	}
	257
	258	}
	259	}
	260
	261	Float_t kZEMcorrPlusPlus[21] = {0.6840,0.7879,0.8722,0.9370,0.9837,1.0137,
	262	1.0292,1.0327,1.0271,1.0152,1.0000,0.9844,
	263	0.9714,0.9634,0.9626,0.9708,0.9891,1.0181,
	264	1.0574,1.1060,1.1617};
	265	Float_t kZEMcorrMoinsMoins[21] = {0.7082,0.8012,0.8809,0.9447,0.9916,1.0220,
	266	1.0372,1.0395,1.0318,1.0174,1.0000,0.9830,
	267	0.9699,0.9635,0.9662,0.9794,1.0034,1.0371,
	268	1.0782,1.1224,1.1634};
	269	///////////////////
	270	//ZemEn correction//
	271	////////////////////
	272
	273	for(Int_t k = -10; k <= 10; ++k)
	274	{
	275	if(zdcAccSat[k+10])
	276	{
	277	//std::cout<<"Vertex at "<<37.5*(Float_t)k<<"cm from cent"<<std::endl;
	278	if(fcurrentDipo>0 && fcurrentL3>0)
	279	{
	280	fzemEn /= kZEMcorrPlusPlus[k+10];
	281	}
	282	if(fcurrentDipo<0 && fcurrentL3<0)
	283	{
	284	fzemEn /= kZEMcorrMoinsMoins[k+10];
	285	}
	286	}
	287	}
	288
	289	////////////////////////
	290	//Centrality selection//
	291	////////////////////////
	292	Double_t fzdcPercentile = -1;
	293	Float_t slope;
294
295
296	if(runnumber < 137722 && runnumber >= 137161 )
297	{
298	if(fzemEn > 345.)
299	{
300	slope = (fzdcEn + 16992.)/(fzemEn - 345.);
301	slope += 2.23902e+02;
302	fzdcPercentile = (TMath::ATan(slope) - 1.56667)/9.49434e-05;
303	if (fzdcPercentile<0) fzdcPercentile = 0;
304	}
305	else fzdcPercentile = 100; }
306	else if(runnumber < 139172 && runnumber >= 137722)
307	{
308	if(fzemEn > 295.)
309	{
310	slope = (fzdcEn + 15000.)/(fzemEn - 295.);
311	slope += 2.23660e+02;
312	fzdcPercentile = (TMath::ATan(slope) - 1.56664)/8.99571e-05;
313	if (fzdcPercentile<0) fzdcPercentile = 0;
314	} else fzdcPercentile = 100; }
315
316
317	else if(runnumber >= 139172)
318	{
319	if(fzemEn > 345.)
320	{
321	slope = (fzdcEn + 16992.)/(fzemEn - 345.);
322	slope += 2.04789e+02;
323	fzdcPercentile = (TMath::ATan(slope) - 1.56629)/1.02768e-04;
324	if (fzdcPercentile<0) fzdcPercentile = 0;
325	} else fzdcPercentile = 100; }
326	else
327	{
328	if(fzemEn > 345.)
329	{
330	slope = (fzdcEn + 16992.)/(fzemEn - 345.);
331	slope += 2.04789e+02;
332	fzdcPercentile = (TMath::ATan(slope) - 1.56629)/1.02768e-04;
333	if(fzdcPercentile<0) fzdcPercentile = 0;
334	}
335	else fzdcPercentile = 100; }
336
337	if(fzdcPercentile > 0 && fzdcPercentile <100) {
338	//std::cout<<fzdcPercentile<<std::endl;
339	cent = fzdcPercentile;
340	}
341	return cent;
342	}
343	//____________________________________________________________________
344	//
345	// EOF
346	//