[u/mrichter/AliRoot.git] / ZDC / AliZDCReconstructor.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

///////////////////////////////////////////////////////////////////////////////
//                                                                           //
// class for ZDC reconstruction                                              //
//                                                                           //
///////////////////////////////////////////////////////////////////////////////


#include <TF1.h>

#include "AliRunLoader.h"
#include "AliRawReader.h"
#include "AliESDEvent.h"
#include "AliZDCDigit.h"
#include "AliZDCRawStream.h"
#include "AliZDCReco.h"
#include "AliZDCReconstructor.h"
#include "AliZDCCalibData.h"


ClassImp(AliZDCReconstructor)


//_____________________________________________________________________________
AliZDCReconstructor:: AliZDCReconstructor() :

  fZNCen (new TF1("fZNCen", 
	"(-2.287920+sqrt(2.287920*2.287920-4*(-0.007629)*(11.921710-x)))/(2*(-0.007629))",0.,164.)),
  fZNPer (new TF1("fZNPer",
      "(-37.812280-sqrt(37.812280*37.812280-4*(-0.190932)*(-1709.249672-x)))/(2*(-0.190932))",0.,164.)),
  fZPCen (new TF1("fZPCen",
       "(-1.321353+sqrt(1.321353*1.321353-4*(-0.007283)*(3.550697-x)))/(2*(-0.007283))",0.,60.)),
  fZPPer (new TF1("fZPPer",
      "(-42.643308-sqrt(42.643308*42.643308-4*(-0.310786)*(-1402.945615-x)))/(2*(-0.310786))",0.,60.)),
  fZDCCen (new TF1("fZDCCen",
      "(-1.934991+sqrt(1.934991*1.934991-4*(-0.004080)*(15.111124-x)))/(2*(-0.004080))",0.,225.)),
  fZDCPer (new TF1("fZDCPer",
      "(-34.380639-sqrt(34.380639*34.380639-4*(-0.104251)*(-2612.189017-x)))/(2*(-0.104251))",0.,225.)),
  fbCen (new TF1("fbCen","-0.056923+0.079703*x-0.0004301*x*x+0.000001366*x*x*x",0.,220.)),
  fbPer (new TF1("fbPer","17.943998-0.046846*x+0.000074*x*x",0.,220.)),
  fZEMn (new TF1("fZEMn","126.2-0.05399*x+0.000005679*x*x",0.,4000.)),
  fZEMp (new TF1("fZEMp","82.49-0.03611*x+0.00000385*x*x",0.,4000.)),
  fZEMsp (new TF1("fZEMsp","208.7-0.09006*x+0.000009526*x*x",0.,4000.)),
  fZEMb (new TF1("fZEMb",
	"16.06-0.01633*x+1.44e-5*x*x-6.778e-9*x*x*x+1.438e-12*x*x*x*x-1.112e-16*x*x*x*x*x",0.,4000.)),
  fCalibData(GetCalibData())

{
  // **** Default constructor

}


//_____________________________________________________________________________
AliZDCReconstructor::~AliZDCReconstructor()
{
// destructor

  delete fZNCen;
  delete fZNPer;
  delete fZPCen;
  delete fZPPer;
  delete fZDCCen;
  delete fZDCPer;
  delete fbCen;
  delete fbPer;
  delete fZEMn;
  delete fZEMp;
  delete fZEMsp;
  delete fZEMb;
}


//_____________________________________________________________________________
void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
{
  // *** Local ZDC reconstruction for digits
  // Works on the current event
    
  Float_t meanPed[47];
  for(Int_t jj=0; jj<47; jj++) meanPed[jj] = fCalibData->GetMeanPed(jj);

  // get digits
  AliZDCDigit digit;
  AliZDCDigit* pdigit = &digit;
  digitsTree->SetBranchAddress("ZDC", &pdigit);

  // loop over digits
  Float_t zn1corr=0, zp1corr=0, zn2corr=0, zp2corr=0, zemcorr=0;
  for (Int_t iDigit = 0; iDigit < digitsTree->GetEntries(); iDigit++) {
    digitsTree->GetEntry(iDigit);
    if (!pdigit) continue;

    if(digit.GetSector(0) == 1)
      zn1corr  += (Float_t) (digit.GetADCValue(0)-meanPed[digit.GetSector(1)]);     // high gain ZN1 ADCs
    else if(digit.GetSector(0) == 2)
      zp1corr  += (Float_t) (digit.GetADCValue(0)-meanPed[digit.GetSector(1)+10]);  // high gain ZP1 ADCs
    else if(digit.GetSector(0) == 3){
      if(digit.GetSector(1)==1)      
	zemcorr += (Float_t) (digit.GetADCValue(0)-meanPed[digit.GetSector(1)+20]); // high gain ZEM1 ADCs
      else if(digit.GetSector(1)==2) 
	zemcorr += (Float_t) (digit.GetADCValue(0)-meanPed[digit.GetSector(1)+22]); // high gain ZEM2 ADCs
    }
    else if(digit.GetSector(0) == 4)
      zn2corr  += (Float_t) (digit.GetADCValue(0)-meanPed[digit.GetSector(1)+24]);  // high gain ZN2 ADCs
    else if(digit.GetSector(0) == 5)
      zp2corr  += (Float_t) (digit.GetADCValue(0)-meanPed[digit.GetSector(1)+34]);  // high gain ZP2 ADCs
  }
  if(zn1corr<0)  zn1corr=0;
  if(zp1corr<0)  zp1corr=0;
  if(zn2corr<0)  zn2corr=0;
  if(zp2corr<0)  zp2corr=0;
  if(zemcorr<0)  zemcorr=0;

  // reconstruct the event
  //printf("\n \t ZDCReco from digits-> ZN = %.0f, ZP = %.0f, ZEM = %.0f\n",zncorr,zpcorr,zemcorr);
  ReconstructEvent(clustersTree, zn1corr, zp1corr, zemcorr, zn2corr, zp2corr);

}

//_____________________________________________________________________________
void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTree) const
{
  // *** ZDC raw data reconstruction
  // Works on the current event
  
  Float_t meanPed[47];
  for(Int_t jj=0; jj<47; jj++) meanPed[jj] = fCalibData->GetMeanPed(jj);

  rawReader->Reset();

  // loop over raw data digits
  Float_t zn1corr=0, zp1corr=0,  zn2corr=0, zp2corr=0,zemcorr=0;
  AliZDCRawStream digit(rawReader);
  while (digit.Next()) {
    if(digit.IsADCDataWord()){
      if(digit.GetADCGain() == 0){
	if(digit.GetSector(0) == 1)	   
	  zn1corr  += (Float_t) (digit.GetADCValue()-meanPed[digit.GetSector(1)]); // high gain ZN1 ADCs;
	else if(digit.GetSector(0) == 2) 
	  zp1corr  += (Float_t) (digit.GetADCValue()-meanPed[digit.GetSector(1)+10]); // high gain ZP1 ADCs;
	else if(digit.GetSector(0) == 3){
	  if(digit.GetSector(1)==1)      
	    zemcorr += (Float_t) (digit.GetADCValue()-meanPed[digit.GetSector(1)+20]); // high gain ZEM1 ADCs
	  else if(digit.GetSector(1)==2) 
	    zemcorr += (Float_t) (digit.GetADCValue()-meanPed[digit.GetSector(1)+22]); // high gain ZEM2 ADCs
	}
	else if(digit.GetSector(0) == 4)	   
	  zn2corr  += (Float_t) (digit.GetADCValue()-meanPed[digit.GetSector(1)+24]); // high gain ZN2 ADCs;
	else if(digit.GetSector(0) == 5) 
	  zp2corr  += (Float_t) (digit.GetADCValue()-meanPed[digit.GetSector(1)+34]); // high gain ZP2 ADCs;
      }
    }
  }
  if(zn1corr<0) zn1corr=0;
  if(zp1corr<0) zp1corr=0;
  if(zn2corr<0) zn2corr=0;
  if(zp2corr<0) zp2corr=0;
  if(zemcorr<0) zemcorr=0;
    
  // reconstruct the event
  //printf("\n\t ZDCReco from raw-> ZN = %.0f, ZP = %.0f, ZEM = %.0f\n",zncorr,zpcorr,zemcorr);
  ReconstructEvent(clustersTree, zn1corr, zp1corr, zemcorr, zn2corr, zp2corr);

}

//_____________________________________________________________________________
void AliZDCReconstructor::ReconstructEvent(TTree *clustersTree, Float_t zn1corr, 
	Float_t zp1corr, Float_t zemcorr, Float_t zn2corr, Float_t zp2corr) const
{
  // ***** Reconstruct one event
  
  //  ---      ADCchannel -> photoelectrons
  // NB-> PM gain = 10^(5), ADC resolution = 6.4*10^(-7)
  // Move to V965 (E.S.,15/09/04) NB-> PM gain = 10^(5), ADC resolution = 8*10^(-7)
  Float_t zn1phe, zp1phe, zemphe, zn2phe, zp2phe, convFactor = 0.08;
  zn1phe  = zn1corr/convFactor;
  zp1phe  = zp1corr/convFactor;
  zemphe = zemcorr/convFactor;
  zn2phe  = zn2corr/convFactor;
  zp2phe  = zp2corr/convFactor;
  //if AliDebug(1,Form("\n    znphe = %f, zpphe = %f, zemphe = %f\n",znphe, zpphe, zemphe);
  
  //  ---      Energy calibration
  // Conversion factors for hadronic ZDCs goes from phe yield to TRUE 
  // incident energy (conversion from GeV to TeV is included); while for EM 
  // calos conversion is from light yield to detected energy calculated by
  // GEANT NB -> ZN and ZP conversion factors are constant since incident
  // spectators have all the same energy, ZEM energy is obtained through a
  // fit over the whole range of incident particle energies 
  // (obtained with full HIJING simulations) 
  Float_t zn1energy, zp1energy, zemenergy, zdc1energy, zn2energy, zp2energy, zdc2energy;
  Float_t zn1phexTeV=329., zp1phexTeV=369., zn2phexTeV=329., zp2phexTeV=369.;
  zn1energy  = zn1phe/zn1phexTeV;
  zp1energy  = zp1phe/zp1phexTeV;
  zdc1energy = zn1energy+zp1energy;
  zn2energy  = zn2phe/zn2phexTeV;
  zp2energy  = zp2phe/zp2phexTeV;
  zdc2energy = zn2energy+zp2energy;
  zemenergy = -4.81+0.3238*zemphe;
  if(zemenergy<0) zemenergy=0;
  //  if AliDebug(1,Form("    znenergy = %f TeV, zpenergy = %f TeV, zdcenergy = %f GeV, "
  //			   "\n		zemenergy = %f TeV\n", znenergy, zpenergy, 
  //			   zdcenergy, zemenergy);
  //  if(zdcenergy==0)
  //    if AliDebug(1,Form("\n\n	###	ATTENZIONE!!! -> ev# %d: znenergy = %f TeV, zpenergy = %f TeV, zdcenergy = %f GeV, "
  //			     " zemenergy = %f TeV\n\n", fMerger->EvNum(), znenergy, zpenergy, zdcenergy, zemenergy); 

  //  ---      Number of detected spectator nucleons
  //  *** N.B. -> It works only in Pb-Pb
  Int_t nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight;
  nDetSpecNLeft = (Int_t) (zn1energy/2.760);
  nDetSpecPLeft = (Int_t) (zp1energy/2.760);
  nDetSpecNRight = (Int_t) (zn2energy/2.760);
  nDetSpecPRight = (Int_t) (zp2energy/2.760);

  //  ---      Number of generated spectator nucleons (from HIJING parameterization)
   //  *** N.B. -> Only one side!!!
 Int_t nGenSpecN=0, nGenSpecP=0, nGenSpec=0;
  Double_t impPar=0;
  // Cut value for Ezem (GeV)
  // ### Results from production  -> 0<b<18 fm (Apr 2002)
  Float_t eZEMCut = 420.;
  Float_t deltaEZEMSup = 690.; 
  Float_t deltaEZEMInf = 270.; 
  if(zemenergy > (eZEMCut+deltaEZEMSup)){
    nGenSpecN = (Int_t) (fZNCen->Eval(zn1energy));
    nGenSpecP = (Int_t) (fZPCen->Eval(zp1energy));
    nGenSpec  = (Int_t) (fZDCCen->Eval(zdc1energy));
    impPar    = fbCen->Eval(zdc1energy);
  }
  else if(zemenergy < (eZEMCut-deltaEZEMInf)){
    nGenSpecN = (Int_t) (fZNPer->Eval(zn1energy)); 
    nGenSpecP = (Int_t) (fZPPer->Eval(zp1energy));
    nGenSpec  = (Int_t) (fZDCPer->Eval(zdc1energy));
    impPar    = fbPer->Eval(zdc1energy);
  }
  else if(zemenergy >= (eZEMCut-deltaEZEMInf) && zemenergy <= (eZEMCut+deltaEZEMSup)){
    nGenSpecN = (Int_t) (fZEMn->Eval(zemenergy));
    nGenSpecP = (Int_t) (fZEMp->Eval(zemenergy));
    nGenSpec  = (Int_t)(fZEMsp->Eval(zemenergy));
    impPar    =  fZEMb->Eval(zemenergy);
  }
  // ### Results from production  -> 0<b<18 fm (Apr 2002)
  if(zn1energy>162.)  nGenSpecN = (Int_t) (fZEMn->Eval(zemenergy));
  if(zp1energy>59.75)  nGenSpecP = (Int_t) (fZEMp->Eval(zemenergy));
  if(zdc1energy>221.5) nGenSpec  = (Int_t)(fZEMsp->Eval(zemenergy));
  if(zdc1energy>220.)  impPar    =  fZEMb->Eval(zemenergy);
  
  if(nGenSpecN>125)    nGenSpecN=125;
  else if(nGenSpecN<0) nGenSpecN=0;
  if(nGenSpecP>82)     nGenSpecP=82;
  else if(nGenSpecP<0) nGenSpecP=0;
  if(nGenSpec>207)     nGenSpec=207;
  else if(nGenSpec<0)  nGenSpec=0;
  
  //  ---      Number of generated participants (from HIJING parameterization)
  Int_t nPart, nPartTot;
  nPart = 207-nGenSpecN-nGenSpecP;
  nPartTot = 207-nGenSpec;
  //printf("\t  ZDCeventReco-> ZNEn = %.0f GeV, ZPEn = %.0f GeV, ZEMEn = %.0f GeV\n",
  //	znenergy, zpenergy, zemenergy);

  // get the output tree and store the ZDC reco object there
  AliZDCReco reco(zn1energy, zp1energy, zdc1energy, zemenergy,
  		  zn2energy, zp2energy, zdc2energy, 
                  nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
		  nGenSpecN, nGenSpecP, nGenSpec,nPartTot, impPar);
  AliZDCReco* preco = &reco;
  const Int_t kBufferSize = 4000;
  clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);

  // write the output tree
  clustersTree->Fill();
}

//_____________________________________________________________________________
void AliZDCReconstructor::FillZDCintoESD(TTree *clustersTree, AliESDEvent* esd) const
{
  // fill energies and number of participants to the ESD

  AliZDCReco reco;
  AliZDCReco* preco = &reco;
  clustersTree->SetBranchAddress("ZDC", &preco);

  clustersTree->GetEntry(0);
  esd->SetZDC(reco.GetZN1energy(), reco.GetZP1energy(), reco.GetZEMenergy(),
	      reco.GetZN2energy(), reco.GetZP2energy(), reco.GetNPart());
}

//_____________________________________________________________________________
AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri) 
{
  // Setting the storage

  Bool_t deleteManager = kFALSE;
  
  AliCDBManager *manager = AliCDBManager::Instance();
  AliCDBStorage *defstorage = manager->GetDefaultStorage();
  
  if(!defstorage || !(defstorage->Contains("ZDC"))){ 
     AliWarning("No default storage set or default storage doesn't contain ZDC!");
     manager->SetDefaultStorage(uri);
     deleteManager = kTRUE;
  }
 
  AliCDBStorage *storage = manager->GetDefaultStorage();

  if(deleteManager){
    AliCDBManager::Instance()->UnsetDefaultStorage();
    defstorage = 0;   // the storage is killed by AliCDBManager::Instance()->Destroy()
  }

  return storage; 
}

//_____________________________________________________________________________
AliZDCCalibData* AliZDCReconstructor::GetCalibData() const
{

  // Getting calibration object for ZDC set

  AliCDBEntry  *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Data");
  if(!entry) AliFatal("No calibration data loaded!");  

  AliZDCCalibData *calibdata = dynamic_cast<AliZDCCalibData*>  (entry->GetObject());
  if(!calibdata)  AliFatal("Wrong calibration object in calibration  file!");

  return calibdata;
}
Commit	Line	Data
8309c1ab	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id$ */
	17
	18	///////////////////////////////////////////////////////////////////////////////
	19	// //
	20	// class for ZDC reconstruction //
	21	// //
	22	///////////////////////////////////////////////////////////////////////////////
	23
	24
	25	#include <TF1.h>
	26
	27	#include "AliRunLoader.h"
	28	#include "AliRawReader.h"
af885e0f	29	#include "AliESDEvent.h"
8309c1ab	30	#include "AliZDCDigit.h"
	31	#include "AliZDCRawStream.h"
	32	#include "AliZDCReco.h"
	33	#include "AliZDCReconstructor.h"
48642b09	34	#include "AliZDCCalibData.h"
8309c1ab	35
	36
	37	ClassImp(AliZDCReconstructor)
	38
	39
	40	//_____________________________________________________________________________
cc2abffd	41	AliZDCReconstructor:: AliZDCReconstructor() :
	42
	43	fZNCen (new TF1("fZNCen",
	44	"(-2.287920+sqrt(2.2879202.287920-4(-0.007629)(11.921710-x)))/(2(-0.007629))",0.,164.)),
	45	fZNPer (new TF1("fZNPer",
	46	"(-37.812280-sqrt(37.81228037.812280-4(-0.190932)(-1709.249672-x)))/(2(-0.190932))",0.,164.)),
	47	fZPCen (new TF1("fZPCen",
	48	"(-1.321353+sqrt(1.3213531.321353-4(-0.007283)(3.550697-x)))/(2(-0.007283))",0.,60.)),
	49	fZPPer (new TF1("fZPPer",
	50	"(-42.643308-sqrt(42.64330842.643308-4(-0.310786)(-1402.945615-x)))/(2(-0.310786))",0.,60.)),
	51	fZDCCen (new TF1("fZDCCen",
	52	"(-1.934991+sqrt(1.9349911.934991-4(-0.004080)(15.111124-x)))/(2(-0.004080))",0.,225.)),
	53	fZDCPer (new TF1("fZDCPer",
	54	"(-34.380639-sqrt(34.38063934.380639-4(-0.104251)(-2612.189017-x)))/(2(-0.104251))",0.,225.)),
	55	fbCen (new TF1("fbCen","-0.056923+0.079703x-0.0004301xx+0.000001366xxx",0.,220.)),
	56	fbPer (new TF1("fbPer","17.943998-0.046846x+0.000074x*x",0.,220.)),
	57	fZEMn (new TF1("fZEMn","126.2-0.05399x+0.000005679x*x",0.,4000.)),
	58	fZEMp (new TF1("fZEMp","82.49-0.03611x+0.00000385x*x",0.,4000.)),
	59	fZEMsp (new TF1("fZEMsp","208.7-0.09006x+0.000009526x*x",0.,4000.)),
	60	fZEMb (new TF1("fZEMb",
	61	"16.06-0.01633x+1.44e-5xx-6.778e-9xxx+1.438e-12xxxx-1.112e-16xxxx*x",0.,4000.)),
	62	fCalibData(GetCalibData())
8309c1ab	63
8309c1ab	64	{
cc2abffd	65	// **** Default constructor
8309c1ab	66
8309c1ab	67	}
8309c1ab	68
8309c1ab	69
	70	//_____________________________________________________________________________
	71	AliZDCReconstructor::~AliZDCReconstructor()
	72	{
	73	// destructor
	74
	75	delete fZNCen;
	76	delete fZNPer;
	77	delete fZPCen;
	78	delete fZPPer;
	79	delete fZDCCen;
	80	delete fZDCPer;
	81	delete fbCen;
	82	delete fbPer;
	83	delete fZEMn;
	84	delete fZEMp;
	85	delete fZEMsp;
	86	delete fZEMb;
	87	}
	88
	89
	90	//_____________________________________________________________________________
70f04f6d	91	void AliZDCReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
8309c1ab	92	{
48642b09	93	// *** Local ZDC reconstruction for digits
70f04f6d	94	// Works on the current event
78d18275	95
48642b09	96	Float_t meanPed[47];
78d18275	97	for(Int_t jj=0; jj<47; jj++) meanPed[jj] = fCalibData->GetMeanPed(jj);
8309c1ab	98
70f04f6d	99	// get digits
8309c1ab	100	AliZDCDigit digit;
8309c1ab	101	AliZDCDigit* pdigit = &digit;
70f04f6d	102	digitsTree->SetBranchAddress("ZDC", &pdigit);
	103
	104	// loop over digits
	105	Float_t zn1corr=0, zp1corr=0, zn2corr=0, zp2corr=0, zemcorr=0;
	106	for (Int_t iDigit = 0; iDigit < digitsTree->GetEntries(); iDigit++) {
	107	digitsTree->GetEntry(iDigit);
	108	if (!pdigit) continue;
	109
	110	if(digit.GetSector(0) == 1)
	111	zn1corr += (Float_t) (digit.GetADCValue(0)-meanPed[digit.GetSector(1)]); // high gain ZN1 ADCs
	112	else if(digit.GetSector(0) == 2)
	113	zp1corr += (Float_t) (digit.GetADCValue(0)-meanPed[digit.GetSector(1)+10]); // high gain ZP1 ADCs
	114	else if(digit.GetSector(0) == 3){
	115	if(digit.GetSector(1)==1)
	116	zemcorr += (Float_t) (digit.GetADCValue(0)-meanPed[digit.GetSector(1)+20]); // high gain ZEM1 ADCs
	117	else if(digit.GetSector(1)==2)
	118	zemcorr += (Float_t) (digit.GetADCValue(0)-meanPed[digit.GetSector(1)+22]); // high gain ZEM2 ADCs
8309c1ab	119	}
70f04f6d	120	else if(digit.GetSector(0) == 4)
	121	zn2corr += (Float_t) (digit.GetADCValue(0)-meanPed[digit.GetSector(1)+24]); // high gain ZN2 ADCs
	122	else if(digit.GetSector(0) == 5)
	123	zp2corr += (Float_t) (digit.GetADCValue(0)-meanPed[digit.GetSector(1)+34]); // high gain ZP2 ADCs
8309c1ab	124	}
70f04f6d	125	if(zn1corr<0) zn1corr=0;
	126	if(zp1corr<0) zp1corr=0;
	127	if(zn2corr<0) zn2corr=0;
	128	if(zp2corr<0) zp2corr=0;
	129	if(zemcorr<0) zemcorr=0;
	130
	131	// reconstruct the event
	132	//printf("\n \t ZDCReco from digits-> ZN = %.0f, ZP = %.0f, ZEM = %.0f\n",zncorr,zpcorr,zemcorr);
	133	ReconstructEvent(clustersTree, zn1corr, zp1corr, zemcorr, zn2corr, zp2corr);
8309c1ab	134
8309c1ab	135	}
	136
	137	//_____________________________________________________________________________
70f04f6d	138	void AliZDCReconstructor::Reconstruct(AliRawReader* rawReader, TTree* clustersTree) const
8309c1ab	139	{
70f04f6d	140	// *** ZDC raw data reconstruction
70f04f6d	141	// Works on the current event
48642b09	142
48642b09	143	Float_t meanPed[47];
78d18275	144	for(Int_t jj=0; jj<47; jj++) meanPed[jj] = fCalibData->GetMeanPed(jj);
8309c1ab	145
70f04f6d	146	rawReader->Reset();
	147
	148	// loop over raw data digits
	149	Float_t zn1corr=0, zp1corr=0, zn2corr=0, zp2corr=0,zemcorr=0;
	150	AliZDCRawStream digit(rawReader);
	151	while (digit.Next()) {
	152	if(digit.IsADCDataWord()){
	153	if(digit.GetADCGain() == 0){
	154	if(digit.GetSector(0) == 1)
	155	zn1corr += (Float_t) (digit.GetADCValue()-meanPed[digit.GetSector(1)]); // high gain ZN1 ADCs;
	156	else if(digit.GetSector(0) == 2)
	157	zp1corr += (Float_t) (digit.GetADCValue()-meanPed[digit.GetSector(1)+10]); // high gain ZP1 ADCs;
	158	else if(digit.GetSector(0) == 3){
	159	if(digit.GetSector(1)==1)
	160	zemcorr += (Float_t) (digit.GetADCValue()-meanPed[digit.GetSector(1)+20]); // high gain ZEM1 ADCs
	161	else if(digit.GetSector(1)==2)
	162	zemcorr += (Float_t) (digit.GetADCValue()-meanPed[digit.GetSector(1)+22]); // high gain ZEM2 ADCs
48642b09	163	}
70f04f6d	164	else if(digit.GetSector(0) == 4)
	165	zn2corr += (Float_t) (digit.GetADCValue()-meanPed[digit.GetSector(1)+24]); // high gain ZN2 ADCs;
	166	else if(digit.GetSector(0) == 5)
	167	zp2corr += (Float_t) (digit.GetADCValue()-meanPed[digit.GetSector(1)+34]); // high gain ZP2 ADCs;
8309c1ab	168	}
8309c1ab	169	}
8309c1ab	170	}
70f04f6d	171	if(zn1corr<0) zn1corr=0;
	172	if(zp1corr<0) zp1corr=0;
	173	if(zn2corr<0) zn2corr=0;
	174	if(zp2corr<0) zp2corr=0;
	175	if(zemcorr<0) zemcorr=0;
	176
	177	// reconstruct the event
	178	//printf("\n\t ZDCReco from raw-> ZN = %.0f, ZP = %.0f, ZEM = %.0f\n",zncorr,zpcorr,zemcorr);
	179	ReconstructEvent(clustersTree, zn1corr, zp1corr, zemcorr, zn2corr, zp2corr);
8309c1ab	180
8309c1ab	181	}
	182
	183	//_____________________________________________________________________________
70f04f6d	184	void AliZDCReconstructor::ReconstructEvent(TTree *clustersTree, Float_t zn1corr,
980685f2	185	Float_t zp1corr, Float_t zemcorr, Float_t zn2corr, Float_t zp2corr) const
8309c1ab	186	{
48642b09	187	// ***** Reconstruct one event
8309c1ab	188
	189	// --- ADCchannel -> photoelectrons
	190	// NB-> PM gain = 10^(5), ADC resolution = 6.4*10^(-7)
	191	// Move to V965 (E.S.,15/09/04) NB-> PM gain = 10^(5), ADC resolution = 8*10^(-7)
980685f2	192	Float_t zn1phe, zp1phe, zemphe, zn2phe, zp2phe, convFactor = 0.08;
	193	zn1phe = zn1corr/convFactor;
	194	zp1phe = zp1corr/convFactor;
8309c1ab	195	zemphe = zemcorr/convFactor;
980685f2	196	zn2phe = zn2corr/convFactor;
980685f2	197	zp2phe = zp2corr/convFactor;
48642b09	198	//if AliDebug(1,Form("\n znphe = %f, zpphe = %f, zemphe = %f\n",znphe, zpphe, zemphe);
8309c1ab	199
	200	// --- Energy calibration
	201	// Conversion factors for hadronic ZDCs goes from phe yield to TRUE
	202	// incident energy (conversion from GeV to TeV is included); while for EM
	203	// calos conversion is from light yield to detected energy calculated by
	204	// GEANT NB -> ZN and ZP conversion factors are constant since incident
	205	// spectators have all the same energy, ZEM energy is obtained through a
	206	// fit over the whole range of incident particle energies
	207	// (obtained with full HIJING simulations)
980685f2	208	Float_t zn1energy, zp1energy, zemenergy, zdc1energy, zn2energy, zp2energy, zdc2energy;
	209	Float_t zn1phexTeV=329., zp1phexTeV=369., zn2phexTeV=329., zp2phexTeV=369.;
	210	zn1energy = zn1phe/zn1phexTeV;
	211	zp1energy = zp1phe/zp1phexTeV;
	212	zdc1energy = zn1energy+zp1energy;
	213	zn2energy = zn2phe/zn2phexTeV;
	214	zp2energy = zp2phe/zp2phexTeV;
	215	zdc2energy = zn2energy+zp2energy;
8309c1ab	216	zemenergy = -4.81+0.3238*zemphe;
8309c1ab	217	if(zemenergy<0) zemenergy=0;
48642b09	218	// if AliDebug(1,Form(" znenergy = %f TeV, zpenergy = %f TeV, zdcenergy = %f GeV, "
	219	// "\n zemenergy = %f TeV\n", znenergy, zpenergy,
	220	// zdcenergy, zemenergy);
48642b09	221	// if(zdcenergy==0)
	222	// if AliDebug(1,Form("\n\n ### ATTENZIONE!!! -> ev# %d: znenergy = %f TeV, zpenergy = %f TeV, zdcenergy = %f GeV, "
	223	// " zemenergy = %f TeV\n\n", fMerger->EvNum(), znenergy, zpenergy, zdcenergy, zemenergy);
8309c1ab	224
980685f2	225	// --- Number of detected spectator nucleons
	226	// *** N.B. -> It works only in Pb-Pb
	227	Int_t nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight;
	228	nDetSpecNLeft = (Int_t) (zn1energy/2.760);
	229	nDetSpecPLeft = (Int_t) (zp1energy/2.760);
	230	nDetSpecNRight = (Int_t) (zn2energy/2.760);
	231	nDetSpecPRight = (Int_t) (zp2energy/2.760);
	232
	233	// --- Number of generated spectator nucleons (from HIJING parameterization)
	234	// *** N.B. -> Only one side!!!
	235	Int_t nGenSpecN=0, nGenSpecP=0, nGenSpec=0;
8309c1ab	236	Double_t impPar=0;
8309c1ab	237	// Cut value for Ezem (GeV)
980685f2	238	// ### Results from production -> 0<b<18 fm (Apr 2002)
8309c1ab	239	Float_t eZEMCut = 420.;
	240	Float_t deltaEZEMSup = 690.;
	241	Float_t deltaEZEMInf = 270.;
	242	if(zemenergy > (eZEMCut+deltaEZEMSup)){
980685f2	243	nGenSpecN = (Int_t) (fZNCen->Eval(zn1energy));
	244	nGenSpecP = (Int_t) (fZPCen->Eval(zp1energy));
	245	nGenSpec = (Int_t) (fZDCCen->Eval(zdc1energy));
	246	impPar = fbCen->Eval(zdc1energy);
8309c1ab	247	}
8309c1ab	248	else if(zemenergy < (eZEMCut-deltaEZEMInf)){
980685f2	249	nGenSpecN = (Int_t) (fZNPer->Eval(zn1energy));
	250	nGenSpecP = (Int_t) (fZPPer->Eval(zp1energy));
	251	nGenSpec = (Int_t) (fZDCPer->Eval(zdc1energy));
	252	impPar = fbPer->Eval(zdc1energy);
8309c1ab	253	}
	254	else if(zemenergy >= (eZEMCut-deltaEZEMInf) && zemenergy <= (eZEMCut+deltaEZEMSup)){
	255	nGenSpecN = (Int_t) (fZEMn->Eval(zemenergy));
	256	nGenSpecP = (Int_t) (fZEMp->Eval(zemenergy));
	257	nGenSpec = (Int_t)(fZEMsp->Eval(zemenergy));
	258	impPar = fZEMb->Eval(zemenergy);
8309c1ab	259	}
980685f2	260	// ### Results from production -> 0<b<18 fm (Apr 2002)
	261	if(zn1energy>162.) nGenSpecN = (Int_t) (fZEMn->Eval(zemenergy));
	262	if(zp1energy>59.75) nGenSpecP = (Int_t) (fZEMp->Eval(zemenergy));
	263	if(zdc1energy>221.5) nGenSpec = (Int_t)(fZEMsp->Eval(zemenergy));
	264	if(zdc1energy>220.) impPar = fZEMb->Eval(zemenergy);
8309c1ab	265
	266	if(nGenSpecN>125) nGenSpecN=125;
	267	else if(nGenSpecN<0) nGenSpecN=0;
	268	if(nGenSpecP>82) nGenSpecP=82;
	269	else if(nGenSpecP<0) nGenSpecP=0;
	270	if(nGenSpec>207) nGenSpec=207;
	271	else if(nGenSpec<0) nGenSpec=0;
8309c1ab	272
980685f2	273	// --- Number of generated participants (from HIJING parameterization)
8309c1ab	274	Int_t nPart, nPartTot;
	275	nPart = 207-nGenSpecN-nGenSpecP;
	276	nPartTot = 207-nGenSpec;
78d18275	277	//printf("\t ZDCeventReco-> ZNEn = %.0f GeV, ZPEn = %.0f GeV, ZEMEn = %.0f GeV\n",
78d18275	278	// znenergy, zpenergy, zemenergy);
8309c1ab	279
70f04f6d	280	// get the output tree and store the ZDC reco object there
980685f2	281	AliZDCReco reco(zn1energy, zp1energy, zdc1energy, zemenergy,
	282	zn2energy, zp2energy, zdc2energy,
	283	nDetSpecNLeft, nDetSpecPLeft, nDetSpecNRight, nDetSpecPRight,
	284	nGenSpecN, nGenSpecP, nGenSpec,nPartTot, impPar);
8309c1ab	285	AliZDCReco* preco = &reco;
8309c1ab	286	const Int_t kBufferSize = 4000;
70f04f6d	287	clustersTree->Branch("ZDC", "AliZDCReco", &preco, kBufferSize);
8309c1ab	288
8309c1ab	289	// write the output tree
70f04f6d	290	clustersTree->Fill();
8309c1ab	291	}
	292
	293	//_____________________________________________________________________________
70f04f6d	294	void AliZDCReconstructor::FillZDCintoESD(TTree clustersTree, AliESDEvent esd) const
8309c1ab	295	{
70f04f6d	296	// fill energies and number of participants to the ESD
8309c1ab	297
8309c1ab	298	AliZDCReco reco;
8309c1ab	299	AliZDCReco* preco = &reco;
70f04f6d	300	clustersTree->SetBranchAddress("ZDC", &preco);
8309c1ab	301
70f04f6d	302	clustersTree->GetEntry(0);
980685f2	303	esd->SetZDC(reco.GetZN1energy(), reco.GetZP1energy(), reco.GetZEMenergy(),
980685f2	304	reco.GetZN2energy(), reco.GetZP2energy(), reco.GetNPart());
8309c1ab	305	}
48642b09	306
48642b09	307	//_____________________________________________________________________________
78d18275	308	AliCDBStorage* AliZDCReconstructor::SetStorage(const char *uri)
48642b09	309	{
cc2abffd	310	// Setting the storage
48642b09	311
78d18275	312	Bool_t deleteManager = kFALSE;
48642b09	313
78d18275	314	AliCDBManager *manager = AliCDBManager::Instance();
78d18275	315	AliCDBStorage *defstorage = manager->GetDefaultStorage();
48642b09	316
78d18275	317	if(!defstorage \|\| !(defstorage->Contains("ZDC"))){
	318	AliWarning("No default storage set or default storage doesn't contain ZDC!");
	319	manager->SetDefaultStorage(uri);
	320	deleteManager = kTRUE;
	321	}
	322
	323	AliCDBStorage *storage = manager->GetDefaultStorage();
	324
	325	if(deleteManager){
	326	AliCDBManager::Instance()->UnsetDefaultStorage();
	327	defstorage = 0; // the storage is killed by AliCDBManager::Instance()->Destroy()
	328	}
	329
	330	return storage;
	331	}
48642b09	332
78d18275	333	//_____________________________________________________________________________
4fda3ba1	334	AliZDCCalibData* AliZDCReconstructor::GetCalibData() const
78d18275	335	{
48642b09	336
4fda3ba1	337	// Getting calibration object for ZDC set
	338
	339	AliCDBEntry *entry = AliCDBManager::Instance()->Get("ZDC/Calib/Data");
457a440d	340	if(!entry) AliFatal("No calibration data loaded!");
4fda3ba1	341
457a440d	342	AliZDCCalibData calibdata = dynamic_cast<AliZDCCalibData> (entry->GetObject());
457a440d	343	if(!calibdata) AliFatal("Wrong calibration object in calibration file!");
48642b09	344
78d18275	345	return calibdata;
48642b09	346	}