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

#include "AliDisplacedVertexSelection.h"
// #include "AliAnalysisManager.h"
#include "AliMCEvent.h"
#include "AliMultiplicity.h"
// #include "AliMCEventHandler.h"
#include "AliHeader.h"
#include "AliGenEventHeader.h"
#include <iostream>
#include <TROOT.h>
#include <TH1D.h>
#include <TH2D.h>
#include "AliESDEvent.h"
#include "AliESDZDC.h"
ClassImp(AliDisplacedVertexSelection)
#if 0
; // This is for Emacs
#endif

//____________________________________________________________________
AliDisplacedVertexSelection::AliDisplacedVertexSelection()
  : TObject(),
    fDeltaTdc(0),
    fSumTdc(0),
    fZdcEnergy(0),
    fZemEnergy(0),
    fCorrelationZemZdc(0),
    fCorrelationSumDelta(0), 
    fVertexZ(kInvalidVtxZ), 
    fCent(100), 
    fHVertexZ(0),
    fHCent(0),
    fMC(kFALSE)
{
}
//____________________________________________________________________
AliDisplacedVertexSelection::AliDisplacedVertexSelection(const AliDisplacedVertexSelection& o)
  : TObject(o),
    fDeltaTdc(o.fDeltaTdc),
    fSumTdc(o.fSumTdc),
    fZdcEnergy(o.fZdcEnergy),
    fZemEnergy(o.fZemEnergy),
    fCorrelationZemZdc(o.fCorrelationZemZdc),
    fCorrelationSumDelta(o.fCorrelationSumDelta), 
    fVertexZ(kInvalidVtxZ), 
    fCent(100), 
    fHVertexZ(0),
    fHCent(0),
    fMC(kFALSE)
{
}
//____________________________________________________________________
AliDisplacedVertexSelection&
AliDisplacedVertexSelection::operator=(const AliDisplacedVertexSelection& o)
{
  if (&o == this) return *this;
  
  fDeltaTdc  = o.fDeltaTdc;
  fSumTdc    = o.fSumTdc;
  fZdcEnergy = o.fZdcEnergy;
  fZemEnergy = o.fZemEnergy;
  fCorrelationZemZdc = o.fCorrelationZemZdc;
  fCorrelationSumDelta = o.fCorrelationSumDelta;
  fMC = o.fMC;
  
  return *this;
}

//____________________________________________________________________
void
AliDisplacedVertexSelection::SetupForData(TList* l, 
					  const char* /* name*/,
					  Bool_t mc)
{
  fMC = mc;
  
  TList* out = new TList;
  out->SetName("displacedVertex");
  out->SetOwner();
  l->Add(out);

  Double_t dVz   = 37.5;
  Double_t vzMin = (-kMaxK-.5) * dVz;
  Double_t vzMax = (+kMaxK+.5) * dVz;

  fHVertexZ = new TH1D("vertexZ", "Interaction point Z", 
		       2*kMaxK+1, vzMin, vzMax);
  fHVertexZ->SetXTitle("IP_{z} [cm]");
  fHVertexZ->SetYTitle("events");
  fHVertexZ->SetDirectory(0);
  fHVertexZ->SetFillColor(kRed+1);
  fHVertexZ->SetFillStyle(3001);
  out->Add(fHVertexZ);

  Int_t    nCent   = 6;
  Double_t bCent[] = { 0, 5, 10, 20, 30, 40, 100 };
  fHCent = new TH1D("cent", "Centrality", nCent, bCent);
  fHCent->SetXTitle("Centrality [%]");
  fHCent->SetYTitle("events");
  fHCent->SetDirectory(0);
  fHCent->SetFillColor(kBlue+1);
  fHCent->SetFillStyle(3001);
  out->Add(fHCent);

  Int_t    nDeltaT   = 2000;
  Double_t minDeltaT = -250;
  Double_t maxDeltaT = +250;
  fDeltaTdc          = new TH1D("DeltaTdc","#DeltaTDC",
				nDeltaT,minDeltaT,maxDeltaT);
  fDeltaTdc->SetXTitle("#DeltaTDC");
  fDeltaTdc->SetDirectory(0);
  out->Add(fDeltaTdc);

  Int_t    nSumT   = nDeltaT;
  Double_t minSumT = minDeltaT;
  Double_t maxSumT = maxDeltaT;
  fSumTdc          = new TH1D("SumTdc","#sumTDC",nSumT,minSumT,maxSumT);
  fSumTdc->SetXTitle("#sumTDC");
  fSumTdc->SetDirectory(0);
  out->Add(fSumTdc);
  
  Int_t    nZdcE   = 10000;
  Double_t minZdcE = 0;
  Double_t maxZdcE = 200000;
  fZdcEnergy           = new TH1D("ZDCEnergy","E_{ZDC}",nZdcE,minZdcE,maxZdcE);
  fZdcEnergy->SetXTitle("E_{ZDC}");
  fZdcEnergy->SetDirectory(0);
  out->Add(fZdcEnergy);
  
  Int_t    nZemE   = 1000;
  Double_t minZemE = -100;
  Double_t maxZemE = 2900;
  fZemEnergy           = new TH1D("ZEMEnergy","E_{ZEM}",nZemE,minZemE,maxZemE);
  fZemEnergy->SetXTitle("E_{ZEM}");
  fZemEnergy->SetDirectory(0);
  out->Add(fZemEnergy);

  fCorrelationZemZdc   = new TH2D("corrZEMZDC","E_{ZEM} vs E_{ZDC}",
				  nZemE, minZemE, maxZemE, 
				  nZdcE, minZdcE, maxZdcE);
  fCorrelationZemZdc->SetXTitle("E_{ZEM}");
  fCorrelationZemZdc->SetYTitle("E_{ZDC}");
  fCorrelationZemZdc->SetDirectory(0);
  out->Add(fCorrelationZemZdc);
  
  fCorrelationSumDelta =  new TH2D("corrSumDelta","#sumTDC vs #DeltaTDC",
				   nSumT, minSumT, maxSumT,
				   nDeltaT, minDeltaT, maxDeltaT);
  fCorrelationSumDelta->SetXTitle("#sum TDC");
  fCorrelationSumDelta->SetYTitle("#DeltaTDC");
  fCorrelationSumDelta->SetDirectory(0);  
  out->Add(fCorrelationSumDelta);

}

  
//____________________________________________________________________
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
}

//____________________________________________________________________
Float_t
AliDisplacedVertexSelection::GetZemCorr(Int_t k, Bool_t minusminus) const
{
  if (-kMaxK > k || k > kMaxK) return 0;

  // Corrections for magnetic fields
  const Float_t kPlusPlus[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 kMoinsMoins[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};
  
  // MC specific code by Hans, is it used - why are ++ and -- the
  // same?
  const Float_t kPlusPlusMC[21]   = { 0.68400, 0.78790, 0.87220,
				      0.93700, 0.98370, 1.08982,
				      1.03518, 1.01534, 1.01840,
				      1.01493, 1.00000, 0.970083,
				      0.979705,0.960576,0.925394,
				      0.92167, 0.971241,1.08338,
				      1.23517, 1.39308, 1.53943 };
  const Float_t kMoinsMoinsMC[21] = { 0.68400, 0.78790, 0.87220,
				      0.9370,  0.98370, 1.08982,
				      1.03518, 1.01534, 1.0184,
				      1.01493, 1.00000, 0.970083,
				      0.979705,0.960576,0.925394,
				      0.92167, 0.971241,1.08338,
				      1.23517, 1.39308, 1.53943 };
  Int_t i = k+kMaxK;
  if (fMC) 
    return minusminus ? kMoinsMoinsMC[i] : kPlusPlusMC[i];
  return minusminus ? kMoinsMoins[i] : kPlusPlus[i];
}


//____________________________________________________________________
Bool_t
AliDisplacedVertexSelection::CheckOutlier(Int_t ivtx, 
					  const AliESDEvent* esd) const
{
  if (fMC) return false;
  if (ivtx <= 4) return false; // Not outliers

  // --- Find outliers -----------------------------------------------
  AliESDVZERO* esdV0 = esd->GetVZEROData(); 
  Double_t nClusters = esd->GetMultiplicity()->GetNumberOfITSClusters(1);

  // parameter from the fit of VZERO multiplicity vs N SPD cluster
  // distribution
  const Double_t slp[8][21] = {{0.000,0.000,0.000,0.000,0.000,0.000,0.000,//0
				0.000,0.000,0.000,0.000,0.000,0.000,0.000,
				0.000,0.000,0.000,0.000,0.000,0.000,0.000},
			       {0.000,0.000,0.000,0.000,0.000,0.000,0.000,//1
				0.000,0.000,0.000,0.000,0.000,0.000,0.000,
				0.000,0.000,0.000,0.000,0.000,0.000,0.000},
			       {0.000,0.000,0.000,0.000,0.000,0.000,0.000,//2
				0.000,0.000,0.000,0.000,0.000,0.000,0.000,
				0.000,0.000,0.000,0.000,0.000,0.000,0.000},
			       {0.000,0.000,0.000,0.000,0.000,0.000,0.000,//3
				0.000,0.000,0.000,0.000,0.000,0.000,0.000,
				0.000,0.000,0.000,0.000,0.000,0.000,0.000},
			       {0.000,0.000,0.000,0.000,0.000,0.155,0.173,//4 
				0.238,0.348,0.384,0.209,0.362,0.433,0.435,
				0.421,0.400,0.403,0.398,0.407,0.394,0.369},
			       {0.000,0.000,0.000,0.000,0.000,0.356,0.402,//5
				0.504,0.650,0.643,0.331,0.543,0.640,0.642,
				0.614,0.591,0.587,0.564,0.572,0.524,0.576},
			       {0.000,0.000,0.000,0.000,0.000,0.386,0.500,//6
				0.692,0.856,0.810,0.390,0.619,0.713,0.708,
				0.670,0.633,0.622,0.593,0.587,0.538,0.587},
			       {0.000,0.000,0.000,0.000,0.000,0.353,0.454,//7
				0.697,0.944,0.941,0.444,0.670,0.746,0.736,
				0.683,0.642,0.619,0.583,0.576,0.535,0.581}};
  const Double_t cst[8][21] = {{0.000,0.000,0.000,0.000,0.000,0.000,0.000,//0
				0.000,0.000,0.000,0.000,0.000,0.000,0.000,
				0.0000,0.000,0.000,0.000,0.000,0.000,0.00},
			       {0.000,0.000,0.000,0.000,0.000,0.000,0.000,//1
				0.000,0.000,0.000,0.000,0.000,0.000,0.000,
				0.0000,0.000,0.000,0.000,0.000,0.000,0.00},
			       {0.000,0.000,0.000,0.000,0.000,0.000,0.000,//2
				0.000,0.000,0.000,0.000,0.000,0.000,0.000,
				0.0000,0.000,0.000,0.000,0.000,0.000,0.00},
			       {0.000,0.000,0.000,0.000,0.000,0.000,0.000,//3
				0.000,0.000,0.000,0.000,0.000,0.000,0.000,
				0.0000,0.000,0.000,0.000,0.000,0.000,0.00},
			       {0.000,0.000,0.000,0.000,0.000,51.67,54.79,//4
				53.51,23.49,24.60,28.47,27.35,24.47,23.93,
				13.162,12.91,5.305,4.007,-14.34,-9.16,-9.43},
			       {0.000,0.000,0.000,0.000,0.000,136.34,88.84,
				84.85,24.31,27.92,35.64,33.15,32.25,23.07,
				13.746,-6.05,-11.95,1.995,-24.91,-22.55,-37.86},
			       {0.000,0.000,0.000,0.000,0.000,380.42,204.22,
				123.89,34.60,29.26,32.49,30.94,19.67,7.760,
				1.1010,-6.01,-15.66,-3.16,-18.34,-17.33,-22.17},
			       {0.000,0.000,0.000,0.000,0.000,128.00,105.04,
				114.05,15.96,21.67,30.11,30.69,27.66,0.363,-0.511,
				-17.67,-22.40,-11.84,-25.65,-24.29,-24.46}};

  Bool_t outlier = kFALSE;
  // loop over VZERO rings
  for (int iring = 4; iring < 8; iring++) { 
    Double_t multRing = 0;
    for (int iCh=0; iCh<8; iCh++) {
      Int_t idx = iCh+iring*8;
      multRing += esdV0->GetMultiplicity(idx);
    }
    
    // Tigher cut for z=0.  Looser cut for suffician statistics and
    // see saturation effects if present.
    Double_t upFac    = (ivtx == 10 ? .35 : .42);
    Double_t downFac  = (ivtx == 10 ? .20 : ivtx < 8 ? .42 : .26);
    Double_t slpUP    = slp[iring][ivtx]*(1.+upFac); //upper cut
    Double_t slpDOWN  = slp[iring][ivtx]*(1.-downFac); //lower cut
    Double_t constant = cst[iring][ivtx];
    //Cut is applied for rings and vertex of interest
    if ((slpDOWN * nClusters + constant) > multRing || 
	(slpUP   * nClusters + constant) < multRing ) { 
	outlier = kTRUE;
    }
  }

  return outlier;
  //--- End of outlier code ----------------------------------------- 
}
//____________________________________________________________________
Bool_t
AliDisplacedVertexSelection::ProcessMC(const AliMCEvent* mcevent)
{
  if (!fMC) return true;
  if (!mcevent) return false;

  AliMCEvent*        event     = const_cast<AliMCEvent*>(mcevent);
  AliHeader*         header    = event->Header();
  AliGenEventHeader* genHeader = header->GenEventHeader();
  TArrayF            vertex;
  genHeader->PrimaryVertex(vertex);
  Double_t zvtx           = vertex.At(2);
  Double_t intTime        = genHeader->InteractionTime();
  const Double_t kVtx[16] = {-187.5, -150., -112.5, -75., -37.5, 0.,
			     37.5,   75.,  112.5, 150., 187.5, 225.,
			     262.5,  300.,  337.5, 375.};
  Int_t vtxClass          = -1;
  Float_t SigmaVtx        = 5.3;
  Float_t SigmaTime       = 0.18;
  
  for(Int_t iVtx=0; iVtx<16; ++iVtx) {
    // Do not do something for nominal
    if (iVtx == 5) continue;
    if ((zvtx >= kVtx[iVtx] - 3.54 * SigmaVtx && 
	 zvtx <= kVtx[iVtx] + 3.54 * SigmaVtx) && 
	(intTime >= (iVtx * 1.25e-9 - 6.25e-9 - 4 * SigmaTime) && 
	 intTime <= (iVtx * 1.25e-9 - 6.25e-9 + 4 * SigmaTime)))  {
      vtxClass = iVtx;
      // break here to not do more 
      break;
    }
  }
  if(vtxClass > -1) {
    fVertexZ = kVtx[vtxClass];
    return true;
  }
  return false;
}
//____________________________________________________________________
Bool_t
AliDisplacedVertexSelection::Process(const AliESDEvent* esd)
{
  fVertexZ = kInvalidVtxZ; // 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;  
  
  // --- 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.;
 
  // HHD/Maxime inclusion - MC check!
  if (fMC) {
    zdcEn      = (2.9 * esdZDC->GetZDCN1Energy() + 
		  7.2 * esdZDC->GetZDCP1Energy() + 
		  3   * esdZDC->GetZDCN2Energy() +
		  8.7 * esdZDC->GetZDCP2Energy());
    zemEn      = 0.57 * (esdZDC->GetZDCEMEnergy(0) + 
			 esdZDC->GetZDCEMEnergy(1));
  }

  // --- 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;
	  }
	}
      }
    }
  }
  fDeltaTdc->Fill(deltaTdc);
  fSumTdc->Fill(sumTdc);
  fCorrelationSumDelta->Fill(sumTdc, deltaTdc);

  // --- Find the vertex ---------------------------------------------
  Int_t ivtx = -1;
  if(deltaTdc!=0. || sumTdc!=0.) {
    Double_t fillVz = kInvalidVtxZ;
    for (Int_t k = -kMaxK; k <= kMaxK; ++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) continue;
      if (k == 0) { 
	fillVz = 0;
	continue;
      }
      // Set the vertex 
      fVertexZ = 37.5 * k;
      fillVz   = fVertexZ;    
      ivtx     = k + 10; // Used for outlier calculations
      // Correct zem energy 
      if(currentDipo>0 && currentL3>0) zemEn /= GetZemCorr(k,false);
      if(currentDipo<0 && currentL3<0) zemEn /= GetZemCorr(k,true);
      // We need to break here, or we could possibly update zemEn again 
      break;
    }
    if (fillVz != kInvalidVtxZ) fHVertexZ->Fill(fillVz);
  }
  fZemEnergy->Fill(zemEn);
  fZdcEnergy->Fill(zdcEn);
  fCorrelationZemZdc->Fill(zemEn, zdcEn);

  // --- Check if this is an outlier event ---------------------------
  if (CheckOutlier(ivtx, esd)) {
    fVertexZ = kInvalidVtxZ;
    return false;
  }

  // --- Calculate the centrality ------------------------------------
  Float_t c1, c2, c3, c4, c5;
  Int_t runnumber = esd->GetRunNumber();
#if 1 // New centrality determination by HHD
  if (runnumber < 137165 ) {
    // ref 137161
    c1 = 16992;
    c2 = 345;
    c3 = 2.23902e+02;
    c4 = 1.56667;
    c5 = 9.49434e-05;
  }
  else if (runnumber <= 137848 && runnumber >= 137230) {
    // ref 137722
    c1  = 15000;
    c2  = 295;
    c3  = 2.23660e+02;
    c4  = 1.56664;
    c5  = 8.99571e-05;
  }
  else if (runnumber <= 138275 && runnumber >= 138125) {
    // ref 137161, yes that's the correct run!
    c1 = 16992;
    c2 = 345;
    c3 = 2.23902e+02;
    c4 = 1.56667;
    c5 = 9.49434e-05;
  }
  else { // if (runnumber <= 139517 && runnumber >= 138358) {
    // Ref run 139172
    c1 = 16992;
    c2 = 345;
    c3 = 2.04789e+02;
    c4 = 1.56629;
    c5 = 1.02768e-04;
  }
#else // Old centrality determination
  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;
  }
#endif
  if (zemEn > c2) { 
    Float_t slope   = (zdcEn + c1) / (zemEn - c2) + c3;
    Float_t zdcCent = (TMath::ATan(slope) - c4) / c5;
    if (zdcCent >= 0) fCent = zdcCent;
  }
  
  fHCent->Fill(fCent);

  return true;
}
//____________________________________________________________________
//
// EOF
//
Commit	Line	Data
65abd48b	1	#include "AliDisplacedVertexSelection.h"
bfab35d9	2	// #include "AliAnalysisManager.h"
	3	#include "AliMCEvent.h"
	4	#include "AliMultiplicity.h"
	5	// #include "AliMCEventHandler.h"
	6	#include "AliHeader.h"
	7	#include "AliGenEventHeader.h"
65abd48b	8	#include <iostream>
65abd48b	9	#include <TROOT.h>
8449e3e0	10	#include <TH1D.h>
bfab35d9	11	#include <TH2D.h>
65abd48b	12	#include "AliESDEvent.h"
	13	#include "AliESDZDC.h"
	14	ClassImp(AliDisplacedVertexSelection)
	15	#if 0
	16	; // This is for Emacs
	17	#endif
	18
	19	//____________________________________________________________________
	20	AliDisplacedVertexSelection::AliDisplacedVertexSelection()
bfab35d9	21	: TObject(),
	22	fDeltaTdc(0),
	23	fSumTdc(0),
	24	fZdcEnergy(0),
	25	fZemEnergy(0),
	26	fCorrelationZemZdc(0),
	27	fCorrelationSumDelta(0),
8449e3e0	28	fVertexZ(kInvalidVtxZ),
	29	fCent(100),
	30	fHVertexZ(0),
bfab35d9	31	fHCent(0),
bfab35d9	32	fMC(kFALSE)
65abd48b	33	{
	34	}
	35	//____________________________________________________________________
	36	AliDisplacedVertexSelection::AliDisplacedVertexSelection(const AliDisplacedVertexSelection& o)
bfab35d9	37	: TObject(o),
	38	fDeltaTdc(o.fDeltaTdc),
	39	fSumTdc(o.fSumTdc),
	40	fZdcEnergy(o.fZdcEnergy),
	41	fZemEnergy(o.fZemEnergy),
	42	fCorrelationZemZdc(o.fCorrelationZemZdc),
	43	fCorrelationSumDelta(o.fCorrelationSumDelta),
8449e3e0	44	fVertexZ(kInvalidVtxZ),
	45	fCent(100),
	46	fHVertexZ(0),
bfab35d9	47	fHCent(0),
bfab35d9	48	fMC(kFALSE)
65abd48b	49	{
	50	}
	51	//____________________________________________________________________
	52	AliDisplacedVertexSelection&
	53	AliDisplacedVertexSelection::operator=(const AliDisplacedVertexSelection& o)
	54	{
bfab35d9	55	if (&o == this) return *this;
	56
	57	fDeltaTdc = o.fDeltaTdc;
	58	fSumTdc = o.fSumTdc;
	59	fZdcEnergy = o.fZdcEnergy;
	60	fZemEnergy = o.fZemEnergy;
	61	fCorrelationZemZdc = o.fCorrelationZemZdc;
	62	fCorrelationSumDelta = o.fCorrelationSumDelta;
	63	fMC = o.fMC;
	64
65abd48b	65	return *this;
	66	}
	67
	68	//____________________________________________________________________
	69	void
8449e3e0	70	AliDisplacedVertexSelection::SetupForData(TList* l,
bfab35d9	71	const char* /* name*/,
bfab35d9	72	Bool_t mc)
65abd48b	73	{
bfab35d9	74	fMC = mc;
bfab35d9	75
8449e3e0	76	TList* out = new TList;
	77	out->SetName("displacedVertex");
	78	out->SetOwner();
	79	l->Add(out);
	80
	81	Double_t dVz = 37.5;
	82	Double_t vzMin = (-kMaxK-.5) * dVz;
	83	Double_t vzMax = (+kMaxK+.5) * dVz;
	84
	85	fHVertexZ = new TH1D("vertexZ", "Interaction point Z",
	86	2*kMaxK+1, vzMin, vzMax);
	87	fHVertexZ->SetXTitle("IP_{z} [cm]");
	88	fHVertexZ->SetYTitle("events");
	89	fHVertexZ->SetDirectory(0);
	90	fHVertexZ->SetFillColor(kRed+1);
	91	fHVertexZ->SetFillStyle(3001);
	92	out->Add(fHVertexZ);
	93
	94	Int_t nCent = 6;
	95	Double_t bCent[] = { 0, 5, 10, 20, 30, 40, 100 };
	96	fHCent = new TH1D("cent", "Centrality", nCent, bCent);
	97	fHCent->SetXTitle("Centrality [%]");
	98	fHCent->SetYTitle("events");
	99	fHCent->SetDirectory(0);
	100	fHCent->SetFillColor(kBlue+1);
	101	fHCent->SetFillStyle(3001);
	102	out->Add(fHCent);
	103
bfab35d9	104	Int_t nDeltaT = 2000;
	105	Double_t minDeltaT = -250;
	106	Double_t maxDeltaT = +250;
	107	fDeltaTdc = new TH1D("DeltaTdc","#DeltaTDC",
	108	nDeltaT,minDeltaT,maxDeltaT);
	109	fDeltaTdc->SetXTitle("#DeltaTDC");
	110	fDeltaTdc->SetDirectory(0);
	111	out->Add(fDeltaTdc);
	112
	113	Int_t nSumT = nDeltaT;
	114	Double_t minSumT = minDeltaT;
a6cd6718	115	Double_t maxSumT = maxDeltaT;
bfab35d9	116	fSumTdc = new TH1D("SumTdc","#sumTDC",nSumT,minSumT,maxSumT);
	117	fSumTdc->SetXTitle("#sumTDC");
	118	fSumTdc->SetDirectory(0);
	119	out->Add(fSumTdc);
	120
	121	Int_t nZdcE = 10000;
	122	Double_t minZdcE = 0;
	123	Double_t maxZdcE = 200000;
	124	fZdcEnergy = new TH1D("ZDCEnergy","E_{ZDC}",nZdcE,minZdcE,maxZdcE);
	125	fZdcEnergy->SetXTitle("E_{ZDC}");
	126	fZdcEnergy->SetDirectory(0);
	127	out->Add(fZdcEnergy);
	128
	129	Int_t nZemE = 1000;
	130	Double_t minZemE = -100;
	131	Double_t maxZemE = 2900;
	132	fZemEnergy = new TH1D("ZEMEnergy","E_{ZEM}",nZemE,minZemE,maxZemE);
	133	fZemEnergy->SetXTitle("E_{ZEM}");
	134	fZemEnergy->SetDirectory(0);
	135	out->Add(fZemEnergy);
	136
	137	fCorrelationZemZdc = new TH2D("corrZEMZDC","E_{ZEM} vs E_{ZDC}",
	138	nZemE, minZemE, maxZemE,
	139	nZdcE, minZdcE, maxZdcE);
	140	fCorrelationZemZdc->SetXTitle("E_{ZEM}");
	141	fCorrelationZemZdc->SetYTitle("E_{ZDC}");
	142	fCorrelationZemZdc->SetDirectory(0);
	143	out->Add(fCorrelationZemZdc);
	144
	145	fCorrelationSumDelta = new TH2D("corrSumDelta","#sumTDC vs #DeltaTDC",
	146	nSumT, minSumT, maxSumT,
	147	nDeltaT, minDeltaT, maxDeltaT);
	148	fCorrelationSumDelta->SetXTitle("#sum TDC");
	149	fCorrelationSumDelta->SetYTitle("#DeltaTDC");
	150	fCorrelationSumDelta->SetDirectory(0);
	151	out->Add(fCorrelationSumDelta);
	152
65abd48b	153	}
bfab35d9	154
65abd48b	155
	156	//____________________________________________________________________
	157	void
	158	AliDisplacedVertexSelection::Print(Option_t*) const
	159	{
6ab100ec	160	#if 0
65abd48b	161	char ind[gROOT->GetDirLevel()+1];
	162	for (Int_t i = 0; i < gROOT->GetDirLevel(); i++) ind[i] = ' ';
	163	ind[gROOT->GetDirLevel()] = '\0';
	164	std::cout << std::boolalpha
	165	<< std::noboolalpha << std::endl;
6ab100ec	166	#endif
65abd48b	167	}
241cca4d	168
bfab35d9	169	//____________________________________________________________________
	170	Float_t
	171	AliDisplacedVertexSelection::GetZemCorr(Int_t k, Bool_t minusminus) const
	172	{
bc31b177	173	if (-kMaxK > k \|\| k > kMaxK) return 0;
bfab35d9	174
	175	// Corrections for magnetic fields
	176	const Float_t kPlusPlus[21] = { 0.6840,0.7879,0.8722,
	177	0.9370,0.9837,1.0137,
	178	1.0292,1.0327,1.0271,
	179	1.0152,1.0000,0.9844,
	180	0.9714,0.9634,0.9626,
	181	0.9708,0.9891,1.0181,
	182	1.0574,1.1060,1.1617};
	183	const Float_t kMoinsMoins[21] = { 0.7082,0.8012,0.8809,
	184	0.9447,0.9916,1.0220,
	185	1.0372,1.0395,1.0318,
	186	1.0174,1.0000,0.9830,
	187	0.9699,0.9635,0.9662,
	188	0.9794,1.0034,1.0371,
	189	1.0782,1.1224,1.1634};
	190
	191	// MC specific code by Hans, is it used - why are ++ and -- the
	192	// same?
	193	const Float_t kPlusPlusMC[21] = { 0.68400, 0.78790, 0.87220,
	194	0.93700, 0.98370, 1.08982,
	195	1.03518, 1.01534, 1.01840,
	196	1.01493, 1.00000, 0.970083,
	197	0.979705,0.960576,0.925394,
	198	0.92167, 0.971241,1.08338,
	199	1.23517, 1.39308, 1.53943 };
	200	const Float_t kMoinsMoinsMC[21] = { 0.68400, 0.78790, 0.87220,
	201	0.9370, 0.98370, 1.08982,
	202	1.03518, 1.01534, 1.0184,
	203	1.01493, 1.00000, 0.970083,
	204	0.979705,0.960576,0.925394,
	205	0.92167, 0.971241,1.08338,
	206	1.23517, 1.39308, 1.53943 };
	207	Int_t i = k+kMaxK;
	208	if (fMC)
	209	return minusminus ? kMoinsMoinsMC[i] : kPlusPlusMC[i];
	210	return minusminus ? kMoinsMoins[i] : kPlusPlus[i];
	211	}
	212
	213
	214	//____________________________________________________________________
	215	Bool_t
	216	AliDisplacedVertexSelection::CheckOutlier(Int_t ivtx,
	217	const AliESDEvent* esd) const
	218	{
	219	if (fMC) return false;
	220	if (ivtx <= 4) return false; // Not outliers
	221
	222	// --- Find outliers -----------------------------------------------
	223	AliESDVZERO* esdV0 = esd->GetVZEROData();
	224	Double_t nClusters = esd->GetMultiplicity()->GetNumberOfITSClusters(1);
	225
	226	// parameter from the fit of VZERO multiplicity vs N SPD cluster
	227	// distribution
	228	const Double_t slp[8][21] = {{0.000,0.000,0.000,0.000,0.000,0.000,0.000,//0
	229	0.000,0.000,0.000,0.000,0.000,0.000,0.000,
	230	0.000,0.000,0.000,0.000,0.000,0.000,0.000},
	231	{0.000,0.000,0.000,0.000,0.000,0.000,0.000,//1
	232	0.000,0.000,0.000,0.000,0.000,0.000,0.000,
	233	0.000,0.000,0.000,0.000,0.000,0.000,0.000},
	234	{0.000,0.000,0.000,0.000,0.000,0.000,0.000,//2
	235	0.000,0.000,0.000,0.000,0.000,0.000,0.000,
	236	0.000,0.000,0.000,0.000,0.000,0.000,0.000},
	237	{0.000,0.000,0.000,0.000,0.000,0.000,0.000,//3
238	0.000,0.000,0.000,0.000,0.000,0.000,0.000,
239	0.000,0.000,0.000,0.000,0.000,0.000,0.000},
240	{0.000,0.000,0.000,0.000,0.000,0.155,0.173,//4
241	0.238,0.348,0.384,0.209,0.362,0.433,0.435,
242	0.421,0.400,0.403,0.398,0.407,0.394,0.369},
243	{0.000,0.000,0.000,0.000,0.000,0.356,0.402,//5
244	0.504,0.650,0.643,0.331,0.543,0.640,0.642,
245	0.614,0.591,0.587,0.564,0.572,0.524,0.576},
246	{0.000,0.000,0.000,0.000,0.000,0.386,0.500,//6
247	0.692,0.856,0.810,0.390,0.619,0.713,0.708,
248	0.670,0.633,0.622,0.593,0.587,0.538,0.587},
249	{0.000,0.000,0.000,0.000,0.000,0.353,0.454,//7
250	0.697,0.944,0.941,0.444,0.670,0.746,0.736,
251	0.683,0.642,0.619,0.583,0.576,0.535,0.581}};
252	const Double_t cst[8][21] = {{0.000,0.000,0.000,0.000,0.000,0.000,0.000,//0
253	0.000,0.000,0.000,0.000,0.000,0.000,0.000,
254	0.0000,0.000,0.000,0.000,0.000,0.000,0.00},
255	{0.000,0.000,0.000,0.000,0.000,0.000,0.000,//1
256	0.000,0.000,0.000,0.000,0.000,0.000,0.000,
257	0.0000,0.000,0.000,0.000,0.000,0.000,0.00},
258	{0.000,0.000,0.000,0.000,0.000,0.000,0.000,//2
259	0.000,0.000,0.000,0.000,0.000,0.000,0.000,
260	0.0000,0.000,0.000,0.000,0.000,0.000,0.00},
261	{0.000,0.000,0.000,0.000,0.000,0.000,0.000,//3
262	0.000,0.000,0.000,0.000,0.000,0.000,0.000,
263	0.0000,0.000,0.000,0.000,0.000,0.000,0.00},
264	{0.000,0.000,0.000,0.000,0.000,51.67,54.79,//4
265	53.51,23.49,24.60,28.47,27.35,24.47,23.93,
266	13.162,12.91,5.305,4.007,-14.34,-9.16,-9.43},
267	{0.000,0.000,0.000,0.000,0.000,136.34,88.84,
268	84.85,24.31,27.92,35.64,33.15,32.25,23.07,
269	13.746,-6.05,-11.95,1.995,-24.91,-22.55,-37.86},
270	{0.000,0.000,0.000,0.000,0.000,380.42,204.22,
271	123.89,34.60,29.26,32.49,30.94,19.67,7.760,
272	1.1010,-6.01,-15.66,-3.16,-18.34,-17.33,-22.17},
273	{0.000,0.000,0.000,0.000,0.000,128.00,105.04,
274	114.05,15.96,21.67,30.11,30.69,27.66,0.363,-0.511,
275	-17.67,-22.40,-11.84,-25.65,-24.29,-24.46}};
276
277	Bool_t outlier = kFALSE;
278	// loop over VZERO rings
279	for (int iring = 4; iring < 8; iring++) {
280	Double_t multRing = 0;
281	for (int iCh=0; iCh<8; iCh++) {
282	Int_t idx = iCh+iring*8;
283	multRing += esdV0->GetMultiplicity(idx);
284	}
285
286	// Tigher cut for z=0. Looser cut for suffician statistics and
287	// see saturation effects if present.
288	Double_t upFac = (ivtx == 10 ? .35 : .42);
289	Double_t downFac = (ivtx == 10 ? .20 : ivtx < 8 ? .42 : .26);
290	Double_t slpUP = slp[iring][ivtx]*(1.+upFac); //upper cut
291	Double_t slpDOWN = slp[iring][ivtx]*(1.-downFac); //lower cut
292	Double_t constant = cst[iring][ivtx];
293	//Cut is applied for rings and vertex of interest
294	if ((slpDOWN * nClusters + constant) > multRing \|\|
295	(slpUP * nClusters + constant) < multRing ) {
296	outlier = kTRUE;
297	}
298	}
299
300	return outlier;
301	//--- End of outlier code -----------------------------------------
302	}
303	//____________________________________________________________________
304	Bool_t
305	AliDisplacedVertexSelection::ProcessMC(const AliMCEvent* mcevent)
306	{
307	if (!fMC) return true;
308	if (!mcevent) return false;
309
310	AliMCEvent* event = const_cast<AliMCEvent*>(mcevent);
311	AliHeader* header = event->Header();
312	AliGenEventHeader* genHeader = header->GenEventHeader();
313	TArrayF vertex;
314	genHeader->PrimaryVertex(vertex);
315	Double_t zvtx = vertex.At(2);
316	Double_t intTime = genHeader->InteractionTime();
317	const Double_t kVtx[16] = {-187.5, -150., -112.5, -75., -37.5, 0.,
318	37.5, 75., 112.5, 150., 187.5, 225.,
319	262.5, 300., 337.5, 375.};
320	Int_t vtxClass = -1;
321	Float_t SigmaVtx = 5.3;
322	Float_t SigmaTime = 0.18;
323
324	for(Int_t iVtx=0; iVtx<16; ++iVtx) {
325	// Do not do something for nominal
326	if (iVtx == 5) continue;
327	if ((zvtx >= kVtx[iVtx] - 3.54 * SigmaVtx &&
328	zvtx <= kVtx[iVtx] + 3.54 * SigmaVtx) &&
329	(intTime >= (iVtx * 1.25e-9 - 6.25e-9 - 4 * SigmaTime) &&
330	intTime <= (iVtx * 1.25e-9 - 6.25e-9 + 4 * SigmaTime))) {
331	vtxClass = iVtx;
332	// break here to not do more
333	break;
334	}
335	}
336	if(vtxClass > -1) {
337	fVertexZ = kVtx[vtxClass];
338	return true;
339	}
340	return false;
341	}
65abd48b	342	//____________________________________________________________________
241cca4d	343	Bool_t
	344	AliDisplacedVertexSelection::Process(const AliESDEvent* esd)
	345	{
8449e3e0	346	fVertexZ = kInvalidVtxZ; // Default vertex value
241cca4d	347	fCent = 100; // Default centrality value
	348
	349	// Some constants
	350	const Float_t kZDCrefSum = -66.5;
	351	const Float_t kZDCrefDelta = -2.10;
	352	const Float_t kZDCsigmaSum = 3.25;
	353	const Float_t kZDCsigmaDelta = 2.25;
	354	const Float_t kZDCsigmaSumSat = 2.50;
	355	const Float_t kZDCsigmaDeltaSat = 1.35;
	356
241cca4d	357	// --- Get the ESD object ------------------------------------------
	358	AliESDZDC* esdZDC = esd->GetESDZDC();
	359	if (!esdZDC) {
	360	AliWarning("No ZDC ESD object!");
	361	return false;
	362	}
	363	Double_t currentL3 = esd->GetCurrentL3();
	364	Double_t currentDipo = esd->GetCurrentDip();
	365
	366	// --- ZDC and ZEM energy signal -----------------------------------
	367	Double_t zdcEn = (esdZDC->GetZDCN1Energy()+
	368	esdZDC->GetZDCP1Energy()+
	369	esdZDC->GetZDCN2Energy()+
	370	esdZDC->GetZDCP2Energy());
	371	Double_t zemEn = (esdZDC->GetZDCEMEnergy(0)+
	372	esdZDC->GetZDCEMEnergy(1))/8.;
bfab35d9	373
	374	// HHD/Maxime inclusion - MC check!
	375	if (fMC) {
	376	zdcEn = (2.9 * esdZDC->GetZDCN1Energy() +
	377	7.2 * esdZDC->GetZDCP1Energy() +
	378	3 * esdZDC->GetZDCN2Energy() +
	379	8.7 * esdZDC->GetZDCP2Energy());
	380	zemEn = 0.57 * (esdZDC->GetZDCEMEnergy(0) +
	381	esdZDC->GetZDCEMEnergy(1));
	382	}
241cca4d	383
	384	// --- Calculate DeltaT and sumT -----------------------------------
	385	Double_t deltaTdc = 0;
	386	Double_t sumTdc = 0;
	387
	388	for(Int_t i = 0; i < 4; ++i) {
	389	if(esdZDC->GetZDCTDCData(10,i) != 0) {
	390	Double_t tdcCnoCorr = 0.025*(esdZDC->GetZDCTDCData(10,i)-
	391	esdZDC->GetZDCTDCData(14,i));
	392	Double_t tdcC = esdZDC->GetZDCTDCCorrected(10,i);
	393	for(Int_t j = 0; j < 4; ++j) {
	394	if(esdZDC->GetZDCTDCData(12,j) != 0) {
	395	Double_t tdcAnoCorr = 0.025*(esdZDC->GetZDCTDCData(12,j)-
	396	esdZDC->GetZDCTDCData(14,j));
	397	Double_t tdcA = esdZDC->GetZDCTDCCorrected(12,j);
	398	if(esdZDC->TestBit(AliESDZDC::kCorrectedTDCFilled)) {
	399	deltaTdc = tdcC-tdcA;
	400	sumTdc = tdcC+tdcA;
	401	}
	402	else {
	403	deltaTdc = tdcCnoCorr-tdcAnoCorr;
	404	sumTdc = tdcCnoCorr+tdcAnoCorr;
	405	}
	406	}
	407	}
	408	}
	409	}
bfab35d9	410	fDeltaTdc->Fill(deltaTdc);
	411	fSumTdc->Fill(sumTdc);
	412	fCorrelationSumDelta->Fill(sumTdc, deltaTdc);
241cca4d	413
241cca4d	414	// --- Find the vertex ---------------------------------------------
bfab35d9	415	Int_t ivtx = -1;
241cca4d	416	if(deltaTdc!=0. \|\| sumTdc!=0.) {
8449e3e0	417	Double_t fillVz = kInvalidVtxZ;
8449e3e0	418	for (Int_t k = -kMaxK; k <= kMaxK; ++k) {
241cca4d	419	Float_t zsat = 2.55F * k;
	420	Float_t delta = (k == 0 ? kZDCsigmaDelta : kZDCsigmaDeltaSat);
	421	Float_t sum = (k == 0 ? kZDCsigmaSum : kZDCsigmaSumSat);
	422	Float_t dT = deltaTdc - kZDCrefDelta - zsat;
	423	Float_t sT = sumTdc - kZDCrefSum - zsat;
	424	Float_t check = dT * dT / delta / delta + sT * sT / sum / sum;
8449e3e0	425	if (check > 1.0) continue;
	426	if (k == 0) {
	427	fillVz = 0;
	428	continue;
	429	}
241cca4d	430	// Set the vertex
241cca4d	431	fVertexZ = 37.5 * k;
bfab35d9	432	fillVz = fVertexZ;
bfab35d9	433	ivtx = k + 10; // Used for outlier calculations
241cca4d	434	// Correct zem energy
bfab35d9	435	if(currentDipo>0 && currentL3>0) zemEn /= GetZemCorr(k,false);
	436	if(currentDipo<0 && currentL3<0) zemEn /= GetZemCorr(k,true);
	437	// We need to break here, or we could possibly update zemEn again
	438	break;
241cca4d	439	}
8449e3e0	440	if (fillVz != kInvalidVtxZ) fHVertexZ->Fill(fillVz);
241cca4d	441	}
bfab35d9	442	fZemEnergy->Fill(zemEn);
	443	fZdcEnergy->Fill(zdcEn);
	444	fCorrelationZemZdc->Fill(zemEn, zdcEn);
	445
	446	// --- Check if this is an outlier event ---------------------------
bc31b177	447	if (CheckOutlier(ivtx, esd)) {
	448	fVertexZ = kInvalidVtxZ;
	449	return false;
	450	}
241cca4d	451
	452	// --- Calculate the centrality ------------------------------------
	453	Float_t c1, c2, c3, c4, c5;
	454	Int_t runnumber = esd->GetRunNumber();
bfab35d9	455	#if 1 // New centrality determination by HHD
	456	if (runnumber < 137165 ) {
	457	// ref 137161
	458	c1 = 16992;
	459	c2 = 345;
	460	c3 = 2.23902e+02;
	461	c4 = 1.56667;
	462	c5 = 9.49434e-05;
	463	}
	464	else if (runnumber <= 137848 && runnumber >= 137230) {
	465	// ref 137722
	466	c1 = 15000;
	467	c2 = 295;
	468	c3 = 2.23660e+02;
	469	c4 = 1.56664;
	470	c5 = 8.99571e-05;
	471	}
	472	else if (runnumber <= 138275 && runnumber >= 138125) {
	473	// ref 137161, yes that's the correct run!
	474	c1 = 16992;
	475	c2 = 345;
	476	c3 = 2.23902e+02;
	477	c4 = 1.56667;
	478	c5 = 9.49434e-05;
	479	}
	480	else { // if (runnumber <= 139517 && runnumber >= 138358) {
	481	// Ref run 139172
	482	c1 = 16992;
	483	c2 = 345;
	484	c3 = 2.04789e+02;
	485	c4 = 1.56629;
	486	c5 = 1.02768e-04;
	487	}
	488	#else // Old centrality determination
241cca4d	489	if (runnumber < 137722 && runnumber >= 137161 ) {
	490	c1 = 16992;
	491	c2 = 345;
	492	c3 = 2.23902e+02;
	493	c4 = 1.56667;
	494	c5 = 9.49434e-05;
	495	}
	496	else if (runnumber < 139172 && runnumber >= 137722) {
	497	c1 = 15000;
	498	c2 = 295;
	499	c3 = 2.23660e+02;
	500	c4 = 1.56664;
	501	c5 = 8.99571e-05;
	502	}
	503	else if (runnumber >= 139172) {
	504	c1 = 16992;
	505	c2 = 345;
	506	c3 = 2.04789e+02;
	507	c4 = 1.56629;
	508	c5 = 1.02768e-04;
	509	}
	510	else {
	511	c1 = 16992;
	512	c2 = 345;
	513	c3 = 2.04789e+02;
	514	c4 = 1.56629;
	515	c5 = 1.02768e-04;
	516	}
bfab35d9	517	#endif
241cca4d	518	if (zemEn > c2) {
	519	Float_t slope = (zdcEn + c1) / (zemEn - c2) + c3;
	520	Float_t zdcCent = (TMath::ATan(slope) - c4) / c5;
	521	if (zdcCent >= 0) fCent = zdcCent;
	522	}
bfab35d9	523
8449e3e0	524	fHCent->Fill(fCent);
65abd48b	525
241cca4d	526	return true;
241cca4d	527	}
65abd48b	528	//____________________________________________________________________
	529	//
	530	// EOF
	531	//