[u/mrichter/AliRoot.git] / PHOS / AliPHOSReconstructor.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$ */

//_________________________________________________________________________
//--
//-- Yves Schutz (SUBATECH) 
// Reconstruction class. Redesigned from the old AliReconstructionner class and 
// derived from STEER/AliReconstructor. 
// 
// --- ROOT system ---

// --- Standard library ---

// --- AliRoot header files ---
#include "AliLog.h"
#include "AliAltroMapping.h"
#include "AliESDEvent.h"
#include "AliESDCaloCluster.h"
#include "AliESDCaloCells.h"
#include "AliPHOSReconstructor.h"
#include "AliPHOSClusterizerv1.h"
#include "AliPHOSTrackSegmentMakerv1.h"
#include "AliPHOSPIDv1.h"
#include "AliPHOSTracker.h"
#include "AliRawReader.h"
#include "AliPHOSTrigger.h"
#include "AliPHOSGeometry.h"
#include "AliPHOSRecoParam.h"
#include "AliPHOSRecoParamEmc.h"
#include "AliPHOSRecoParamCpv.h"
#include "AliPHOSDigit.h"
#include "AliPHOSTrackSegment.h"
#include "AliPHOSEmcRecPoint.h"
#include "AliPHOSRecParticle.h"
#include "AliPHOSRawDecoder.h"
#include "AliPHOSRawDecoderv1.h"
#include "AliPHOSRawDecoderv2.h"
#include "AliPHOSRawDigiProducer.h"
#include "AliPHOSPulseGenerator.h"

ClassImp(AliPHOSReconstructor)

Bool_t AliPHOSReconstructor::fgDebug = kFALSE ; 
AliPHOSRecoParam* AliPHOSReconstructor::fgkRecoParamEmc =0;  // EMC rec. parameters
AliPHOSRecoParam* AliPHOSReconstructor::fgkRecoParamCpv =0;  // CPV rec. parameters

//____________________________________________________________________________
AliPHOSReconstructor::AliPHOSReconstructor() :
  fGeom(NULL)
{
  // ctor

  if (!fgkRecoParamEmc) {
    AliWarning("The Reconstruction parameters for EMC nonitialized - Used default one");
    fgkRecoParamEmc = AliPHOSRecoParamEmc::GetEmcDefaultParameters();
  }

  if (!fgkRecoParamCpv) {
    AliWarning("The Reconstruction parameters for CPV nonitialized - Used default one");
    fgkRecoParamCpv = AliPHOSRecoParamCpv::GetCpvDefaultParameters();
  }

  fGeom = AliPHOSGeometry::GetInstance("IHEP","");
}

//____________________________________________________________________________
  AliPHOSReconstructor::~AliPHOSReconstructor()
{
  // dtor
  delete fGeom;
} 

//____________________________________________________________________________
void AliPHOSReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
{
  // 'single-event' local reco method called by AliReconstruction; 
  // Only the clusterization is performed,; the rest of the reconstruction is done in FillESD because the track
  // segment maker needs access to the AliESDEvent object to retrieve the tracks reconstructed by 
  // the global tracking.

  AliPHOSClusterizerv1 clu(fGeom);
  clu.SetInput(digitsTree);
  clu.SetOutput(clustersTree);
  if ( Debug() ) 
    clu.Digits2Clusters("deb all") ; 
  else 
    clu.Digits2Clusters("") ;
}

//____________________________________________________________________________
void AliPHOSReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree, 
				   AliESDEvent* esd) const
{
  // This method produces PHOS rec-particles,
  // then it creates AliESDtracks out of them and
  // write tracks to the ESD

  AliPHOSTrackSegmentMaker *tsm = new AliPHOSTrackSegmentMakerv1(fGeom);
  AliPHOSPID               *pid = new AliPHOSPIDv1              (fGeom);

  // do current event; the loop over events is done by AliReconstruction::Run()
  tsm->SetESD(esd) ; 
  tsm->SetInput(clustersTree);
  if ( Debug() ) 
    tsm->Clusters2TrackSegments("deb all") ;
  else 
    tsm->Clusters2TrackSegments("") ;
  
  pid->SetInput(clustersTree, tsm->GetTrackSegments()) ; 
  pid->SetESD(esd) ; 
  if ( Debug() ) 
    pid->TrackSegments2RecParticles("deb all") ;
  else 
    pid->TrackSegments2RecParticles("") ;


  // This function creates AliESDtracks from AliPHOSRecParticles
  //         and
  // writes them to the ESD

  TClonesArray *recParticles  = pid->GetRecParticles();
  Int_t nOfRecParticles = recParticles->GetEntries();
  
  esd->SetNumberOfPHOSClusters(nOfRecParticles) ; 
  esd->SetFirstPHOSCluster(esd->GetNumberOfCaloClusters()) ;
  
  AliDebug(2,Form("%d rec. particles, option %s",nOfRecParticles,GetOption()));
  
  // Read digits array

  TBranch *branch = digitsTree->GetBranch("PHOS");
  if (!branch) { 
    AliError("can't get the branch with the PHOS digits !");
    return;
  }
  TClonesArray *digitsArray    = new TClonesArray("AliPHOSDigit",100);
  branch->SetAddress(&digitsArray);
  branch->GetEntry(0);

  // Get the clusters array

  TBranch *emcbranch = clustersTree->GetBranch("PHOSEmcRP");
  if (!emcbranch) { 
    AliError("can't get the branch with the PHOS EMC clusters !");
    return;
  }

  TObjArray *emcRecPoints = new TObjArray(100) ;
  emcbranch->SetAddress(&emcRecPoints);
  emcbranch->GetEntry(0);

  //#########Calculate trigger and set trigger info###########

  AliPHOSTrigger tr ;
  //   tr.SetPatchSize(1);//create 4x4 patches
  tr.SetSimulation(kFALSE);
  tr.Trigger(digitsArray);
  
  Float_t maxAmp2x2  = tr.Get2x2MaxAmplitude();
  Float_t maxAmpnxn  = tr.GetnxnMaxAmplitude();
  Float_t ampOutOfPatch2x2  = tr.Get2x2AmpOutOfPatch() ;
  Float_t ampOutOfPatchnxn  = tr.GetnxnAmpOutOfPatch() ;

  Int_t iSM2x2      = tr.Get2x2SuperModule();
  Int_t iSMnxn      = tr.GetnxnSuperModule();
  Int_t iCrystalPhi2x2 = tr.Get2x2CrystalPhi();
  Int_t iCrystalPhinxn = tr.GetnxnCrystalPhi();
  Int_t iCrystalEta2x2 = tr.Get2x2CrystalEta();
  Int_t iCrystalEtanxn = tr.GetnxnCrystalEta();

  AliDebug(2, Form("Trigger 2x2 max amp %f, out amp %f, SM %d, iphi %d ieta %d",  
		   maxAmp2x2, ampOutOfPatch2x2, iSM2x2,iCrystalPhi2x2, iCrystalEta2x2));
  AliDebug(2, Form("Trigger 4x4 max amp %f , out amp %f, SM %d, iphi %d, ieta %d",
		   maxAmpnxn, ampOutOfPatchnxn, iSMnxn,iCrystalPhinxn, iCrystalEtanxn));

  // Attention! PHOS modules in order to calculate AbsId need to be 1-5 not 0-4 as returns trigger.
  Int_t iRelId2x2 []= {iSM2x2+1,0,iCrystalPhi2x2,iCrystalEta2x2};
  Int_t iAbsId2x2 =-1;
  Int_t iRelIdnxn []= {iSMnxn+1,0,iCrystalPhinxn,iCrystalEtanxn};
  Int_t iAbsIdnxn =-1;
  TVector3    pos2x2(-1,-1,-1);
  TVector3    posnxn(-1,-1,-1);
  fGeom->RelToAbsNumbering(iRelId2x2, iAbsId2x2);
  fGeom->RelToAbsNumbering(iRelIdnxn, iAbsIdnxn);
  fGeom->RelPosInAlice(iAbsId2x2, pos2x2);
  fGeom->RelPosInAlice(iAbsIdnxn, posnxn);

  TArrayF triggerPosition(6);
  triggerPosition[0] = pos2x2(0) ;   
  triggerPosition[1] = pos2x2(1) ;   
  triggerPosition[2] = pos2x2(2) ;  
  triggerPosition[3] = posnxn(0) ;   
  triggerPosition[4] = posnxn(1) ;   
  triggerPosition[5] = posnxn(2) ;  

  TArrayF triggerAmplitudes(4);
  triggerAmplitudes[0] = maxAmp2x2 ;   
  triggerAmplitudes[1] = ampOutOfPatch2x2 ;    
  triggerAmplitudes[2] = maxAmpnxn ;   
  triggerAmplitudes[3] = ampOutOfPatchnxn ;   

  //esd->SetPHOSTriggerCells(triggerPosition);
  esd->AddPHOSTriggerPosition(triggerPosition);
  esd->AddPHOSTriggerAmplitudes(triggerAmplitudes);
  

  //########################################
  //############# Fill CaloCells ###########
  //########################################

  Int_t nDigits = digitsArray->GetEntries();
  Int_t idignew = 0 ;
  AliDebug(1,Form("%d digits",nDigits));

  const Int_t knEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
  AliESDCaloCells &phsCells = *(esd->GetPHOSCells());
  phsCells.CreateContainer(nDigits);
  phsCells.SetType(AliESDCaloCells::kPHOSCell);

  // Add to CaloCells only EMC digits with non-zero energy 
  for (Int_t idig = 0 ; idig < nDigits ; idig++) {
    const AliPHOSDigit * dig = (const AliPHOSDigit*)digitsArray->At(idig);
    if(dig->GetId() <= knEMC && dig->GetEnergy() > 0 ){
      //printf("i %d; id %d; amp %f; time %e\n",
      //idignew,dig->GetId(),dig->GetEnergy(), dig->GetTime());
      phsCells.SetCell(idignew,dig->GetId(), dig->GetEnergy(), dig->GetTime());   
      idignew++;
    }
  }
  phsCells.SetNumberOfCells(idignew);
  phsCells.Sort();

  //########################################
  //############## Fill CaloClusters #######
  //########################################

  for (Int_t recpart = 0 ; recpart < nOfRecParticles ; recpart++) {
    AliPHOSRecParticle  *rp    = dynamic_cast<AliPHOSRecParticle*>(recParticles->At(recpart));
    if (Debug()) 
      rp->Print();
    // Get track segment and EMC rec.point associated with this rec.particle
    AliPHOSTrackSegment *ts    = static_cast<AliPHOSTrackSegment *>(tsm->GetTrackSegments()
								    ->At(rp->GetPHOSTSIndex()));

    AliPHOSEmcRecPoint  *emcRP = static_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(ts->GetEmcIndex()));
    AliESDCaloCluster   *ec    = new AliESDCaloCluster() ; 
    
    Float_t xyz[3];
    for (Int_t ixyz=0; ixyz<3; ixyz++) 
      xyz[ixyz] = rp->GetPos()[ixyz];
    
    AliDebug(2,Form("Global position xyz=(%f,%f,%f)",xyz[0],xyz[1],xyz[2]));
   
    // Create cell lists

    Int_t     cellMult   = emcRP->GetDigitsMultiplicity();
    Int_t    *digitsList = emcRP->GetDigitsList();
    Float_t  *rpElist    = emcRP->GetEnergiesList() ;
    UShort_t *absIdList  = new UShort_t[cellMult];
    Double_t *fracList   = new Double_t[cellMult];

    for (Int_t iCell=0; iCell<cellMult; iCell++) {
      AliPHOSDigit *digit = static_cast<AliPHOSDigit *>(digitsArray->At(digitsList[iCell]));
      absIdList[iCell] = (UShort_t)(digit->GetId());
      if (digit->GetEnergy() > 0)
 	fracList[iCell] = rpElist[iCell]/digit->GetEnergy();
      else
 	fracList[iCell] = 0;
    }

    //Primaries
    Int_t  primMult  = 0;
    Int_t *primList =  emcRP->GetPrimaries(primMult);

    // fills the ESDCaloCluster
    ec->SetClusterType(AliESDCaloCluster::kPHOSCluster);
    ec->SetPosition(xyz);                       //rec.point position in MARS
    ec->SetE(rp->Energy());                     //total particle energy
    ec->SetClusterDisp(emcRP->GetDispersion()); //cluster dispersion
    ec->SetPid(rp->GetPID()) ;                  //array of particle identification
    ec->SetM02(emcRP->GetM2x()) ;               //second moment M2x
    ec->SetM20(emcRP->GetM2z()) ;               //second moment M2z
    ec->SetNExMax(emcRP->GetNExMax());          //number of local maxima
    ec->SetEmcCpvDistance(ts->GetCpvDistance("r")); //Only radius, what about separate x,z????
    ec->SetClusterChi2(-1);                     //not yet implemented
    ec->SetM11(-1) ;                            //not yet implemented
 
    //Cells contributing to clusters
    ec->SetNCells(cellMult);
    ec->SetCellsAbsId(absIdList);
    ec->SetCellsAmplitudeFraction(fracList);

    //Distance to the nearest bad crystal
    ec->SetDistanceToBadChannel(emcRP->GetDistanceToBadCrystal()); 
  
    //Array of MC indeces
    TArrayI arrayPrim(primMult,primList);
    ec->AddLabels(arrayPrim);

    //Array of tracks uncomment when available in future
    //TArrayS arrayTrackMatched(1);// Only one track, temporal solution.
    //arrayTrackMatched[0]= (Short_t)(matchedTrack[iClust]);
    //ec->AddTracksMatched(arrayTrackMatched);
    
    // add the track to the esd object

    esd->AddCaloCluster(ec);
    delete ec;   
  }
  digitsArray ->Delete();
  delete digitsArray;
  emcRecPoints->Delete();
  delete emcRecPoints;
  delete tsm;
  delete pid;
}

//____________________________________________________________________________
AliTracker* AliPHOSReconstructor::CreateTracker() const
{
  // creates the PHOS tracker
  return new AliPHOSTracker();
}

//____________________________________________________________________________
void  AliPHOSReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digitsTree) const
{
  // Converts raw data to
  // PHOS digits
  // Works on a single-event basis

  rawReader->Reset() ; 

  AliPHOSRawDecoder * dc ;

  const TObjArray* maps = AliPHOSRecoParamEmc::GetMappings();
  if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");

  AliAltroMapping *mapping[4];
  for(Int_t i = 0; i < 4; i++) {
    mapping[i] = (AliAltroMapping*)maps->At(i);
  }

  if(strcmp(fgkRecoParamEmc->DecoderVersion(),"v1")==0) 
    dc=new AliPHOSRawDecoderv1(rawReader,mapping);
  else 
    if(strcmp(fgkRecoParamEmc->DecoderVersion(),"v2")==0) 
      dc=new AliPHOSRawDecoderv2(rawReader,mapping);
    else
      dc=new AliPHOSRawDecoder(rawReader,mapping);

  dc->SetOldRCUFormat(fgkRecoParamEmc->IsOldRCUFormat());
  
  dc->SubtractPedestals(fgkRecoParamEmc->SubtractPedestals());
  
  TClonesArray *digits = new TClonesArray("AliPHOSDigit",1);
  digits->SetName("DIGITS");
  Int_t bufsize = 32000;
  digitsTree->Branch("PHOS", &digits, bufsize);

  AliPHOSRawDigiProducer pr(fgkRecoParamEmc,fgkRecoParamCpv);
  pr.MakeDigits(digits,dc);

  delete dc ;

  //!!!!for debug!!!
/*
  Int_t modMax=-111;
  Int_t colMax=-111;
  Int_t rowMax=-111;
  Float_t eMax=-333;
  //!!!for debug!!!

  Int_t relId[4];
  for(Int_t iDigit=0; iDigit<digits->GetEntries(); iDigit++) {
    AliPHOSDigit* digit = (AliPHOSDigit*)digits->At(iDigit);
    if(digit->GetEnergy()>eMax) {
      fGeom->AbsToRelNumbering(digit->GetId(),relId);
      eMax=digit->GetEnergy();
      modMax=relId[0];
      rowMax=relId[2];
      colMax=relId[3];
    }
  }

  AliDebug(1,Form("Digit with max. energy:  modMax %d colMax %d rowMax %d  eMax %f\n\n",
		  modMax,colMax,rowMax,eMax));
*/
  digitsTree->Fill();
  digits->Delete();
  delete digits;
}
Commit	Line	Data
d15a28e7	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
b2a60966	16	/* $Id$ */
b2a60966	17
d15a28e7	18	//_________________________________________________________________________
a5e00016	19	//--
a5e00016	20	//-- Yves Schutz (SUBATECH)
dfe0be07	21	// Reconstruction class. Redesigned from the old AliReconstructionner class and
	22	// derived from STEER/AliReconstructor.
	23	//
d15a28e7	24	// --- ROOT system ---
d15a28e7	25
d15a28e7	26	// --- Standard library ---
364de5c6	27
d15a28e7	28	// --- AliRoot header files ---
9a2cdbdf	29	#include "AliLog.h"
b3006690	30	#include "AliAltroMapping.h"
af885e0f	31	#include "AliESDEvent.h"
33ba95c3	32	#include "AliESDCaloCluster.h"
a5e00016	33	#include "AliESDCaloCells.h"
f444a19f	34	#include "AliPHOSReconstructor.h"
7acf6008	35	#include "AliPHOSClusterizerv1.h"
7acf6008	36	#include "AliPHOSTrackSegmentMakerv1.h"
7acf6008	37	#include "AliPHOSPIDv1.h"
23904d16	38	#include "AliPHOSTracker.h"
d22dd3b4	39	#include "AliRawReader.h"
64df000d	40	#include "AliPHOSTrigger.h"
64df000d	41	#include "AliPHOSGeometry.h"
86121c91	42	#include "AliPHOSRecoParam.h"
3799bcb5	43	#include "AliPHOSRecoParamEmc.h"
3799bcb5	44	#include "AliPHOSRecoParamCpv.h"
9a2cdbdf	45	#include "AliPHOSDigit.h"
	46	#include "AliPHOSTrackSegment.h"
	47	#include "AliPHOSEmcRecPoint.h"
	48	#include "AliPHOSRecParticle.h"
	49	#include "AliPHOSRawDecoder.h"
77ea1c6f	50	#include "AliPHOSRawDecoderv1.h"
710634f4	51	#include "AliPHOSRawDecoderv2.h"
9a2cdbdf	52	#include "AliPHOSRawDigiProducer.h"
9a2cdbdf	53	#include "AliPHOSPulseGenerator.h"
e957fea8	54
f444a19f	55	ClassImp(AliPHOSReconstructor)
d15a28e7	56
2e60107f	57	Bool_t AliPHOSReconstructor::fgDebug = kFALSE ;
3799bcb5	58	AliPHOSRecoParam* AliPHOSReconstructor::fgkRecoParamEmc =0; // EMC rec. parameters
3799bcb5	59	AliPHOSRecoParam* AliPHOSReconstructor::fgkRecoParamCpv =0; // CPV rec. parameters
2e60107f	60
d15a28e7	61	//____________________________________________________________________________
9a2cdbdf	62	AliPHOSReconstructor::AliPHOSReconstructor() :
9a2cdbdf	63	fGeom(NULL)
d15a28e7	64	{
b2a60966	65	// ctor
0379a13e	66
3799bcb5	67	if (!fgkRecoParamEmc) {
	68	AliWarning("The Reconstruction parameters for EMC nonitialized - Used default one");
	69	fgkRecoParamEmc = AliPHOSRecoParamEmc::GetEmcDefaultParameters();
	70	}
	71
	72	if (!fgkRecoParamCpv) {
	73	AliWarning("The Reconstruction parameters for CPV nonitialized - Used default one");
	74	fgkRecoParamCpv = AliPHOSRecoParamCpv::GetCpvDefaultParameters();
	75	}
	76
9a2cdbdf	77	fGeom = AliPHOSGeometry::GetInstance("IHEP","");
e68222ce	78	}
6ad0bfa0	79
0379a13e	80	//____________________________________________________________________________
	81	AliPHOSReconstructor::~AliPHOSReconstructor()
	82	{
	83	// dtor
9a2cdbdf	84	delete fGeom;
0379a13e	85	}
7acf6008	86
7acf6008	87	//____________________________________________________________________________
9a2cdbdf	88	void AliPHOSReconstructor::Reconstruct(TTree* digitsTree, TTree* clustersTree) const
35293055	89	{
9a2cdbdf	90	// 'single-event' local reco method called by AliReconstruction;
dfe0be07	91	// Only the clusterization is performed,; the rest of the reconstruction is done in FillESD because the track
af885e0f	92	// segment maker needs access to the AliESDEvent object to retrieve the tracks reconstructed by
dfe0be07	93	// the global tracking.
9a2cdbdf	94
	95	AliPHOSClusterizerv1 clu(fGeom);
	96	clu.SetInput(digitsTree);
	97	clu.SetOutput(clustersTree);
a68156e6	98	if ( Debug() )
9a2cdbdf	99	clu.Digits2Clusters("deb all") ;
a68156e6	100	else
9a2cdbdf	101	clu.Digits2Clusters("") ;
a68156e6	102	}
	103
	104	//____________________________________________________________________________
9a2cdbdf	105	void AliPHOSReconstructor::FillESD(TTree* digitsTree, TTree* clustersTree,
9a2cdbdf	106	AliESDEvent* esd) const
a68156e6	107	{
9a2cdbdf	108	// This method produces PHOS rec-particles,
	109	// then it creates AliESDtracks out of them and
	110	// write tracks to the ESD
	111
	112	AliPHOSTrackSegmentMaker *tsm = new AliPHOSTrackSegmentMakerv1(fGeom);
e68222ce	113	AliPHOSPID *pid = new AliPHOSPIDv1 (fGeom);
e68222ce	114
9a2cdbdf	115	// do current event; the loop over events is done by AliReconstruction::Run()
	116	tsm->SetESD(esd) ;
	117	tsm->SetInput(clustersTree);
	118	if ( Debug() )
	119	tsm->Clusters2TrackSegments("deb all") ;
	120	else
	121	tsm->Clusters2TrackSegments("") ;
35293055	122
9a2cdbdf	123	pid->SetInput(clustersTree, tsm->GetTrackSegments()) ;
9a2cdbdf	124	pid->SetESD(esd) ;
dfe0be07	125	if ( Debug() )
9a2cdbdf	126	pid->TrackSegments2RecParticles("deb all") ;
dfe0be07	127	else
9a2cdbdf	128	pid->TrackSegments2RecParticles("") ;
7acf6008	129
77ea1c6f	130
23904d16	131	// This function creates AliESDtracks from AliPHOSRecParticles
	132	// and
	133	// writes them to the ESD
bf72996e	134
9a2cdbdf	135	TClonesArray *recParticles = pid->GetRecParticles();
dfe0be07	136	Int_t nOfRecParticles = recParticles->GetEntries();
a5fa6165	137
85c60a8e	138	esd->SetNumberOfPHOSClusters(nOfRecParticles) ;
dd7ee508	139	esd->SetFirstPHOSCluster(esd->GetNumberOfCaloClusters()) ;
a5e00016	140
9a2cdbdf	141	AliDebug(2,Form("%d rec. particles, option %s",nOfRecParticles,GetOption()));
a5e00016	142
	143	// Read digits array
	144
	145	TBranch *branch = digitsTree->GetBranch("PHOS");
	146	if (!branch) {
	147	AliError("can't get the branch with the PHOS digits !");
	148	return;
	149	}
	150	TClonesArray *digitsArray = new TClonesArray("AliPHOSDigit",100);
	151	branch->SetAddress(&digitsArray);
	152	branch->GetEntry(0);
	153
	154	// Get the clusters array
	155
	156	TBranch *emcbranch = clustersTree->GetBranch("PHOSEmcRP");
	157	if (!emcbranch) {
	158	AliError("can't get the branch with the PHOS EMC clusters !");
	159	return;
	160	}
	161
	162	TObjArray *emcRecPoints = new TObjArray(100) ;
	163	emcbranch->SetAddress(&emcRecPoints);
	164	emcbranch->GetEntry(0);
8013c27e	165
64df000d	166	//#########Calculate trigger and set trigger info###########
a5e00016	167
64df000d	168	AliPHOSTrigger tr ;
64df000d	169	// tr.SetPatchSize(1);//create 4x4 patches
bfae5a5d	170	tr.SetSimulation(kFALSE);
bfae5a5d	171	tr.Trigger(digitsArray);
64df000d	172
	173	Float_t maxAmp2x2 = tr.Get2x2MaxAmplitude();
	174	Float_t maxAmpnxn = tr.GetnxnMaxAmplitude();
	175	Float_t ampOutOfPatch2x2 = tr.Get2x2AmpOutOfPatch() ;
	176	Float_t ampOutOfPatchnxn = tr.GetnxnAmpOutOfPatch() ;
	177
64df000d	178	Int_t iSM2x2 = tr.Get2x2SuperModule();
	179	Int_t iSMnxn = tr.GetnxnSuperModule();
	180	Int_t iCrystalPhi2x2 = tr.Get2x2CrystalPhi();
	181	Int_t iCrystalPhinxn = tr.GetnxnCrystalPhi();
	182	Int_t iCrystalEta2x2 = tr.Get2x2CrystalEta();
	183	Int_t iCrystalEtanxn = tr.GetnxnCrystalEta();
	184
a5e00016	185	AliDebug(2, Form("Trigger 2x2 max amp %f, out amp %f, SM %d, iphi %d ieta %d",
	186	maxAmp2x2, ampOutOfPatch2x2, iSM2x2,iCrystalPhi2x2, iCrystalEta2x2));
	187	AliDebug(2, Form("Trigger 4x4 max amp %f , out amp %f, SM %d, iphi %d, ieta %d",
	188	maxAmpnxn, ampOutOfPatchnxn, iSMnxn,iCrystalPhinxn, iCrystalEtanxn));
64df000d	189
bfae5a5d	190	// Attention! PHOS modules in order to calculate AbsId need to be 1-5 not 0-4 as returns trigger.
bfae5a5d	191	Int_t iRelId2x2 []= {iSM2x2+1,0,iCrystalPhi2x2,iCrystalEta2x2};
64df000d	192	Int_t iAbsId2x2 =-1;
bfae5a5d	193	Int_t iRelIdnxn []= {iSMnxn+1,0,iCrystalPhinxn,iCrystalEtanxn};
64df000d	194	Int_t iAbsIdnxn =-1;
	195	TVector3 pos2x2(-1,-1,-1);
	196	TVector3 posnxn(-1,-1,-1);
9a2cdbdf	197	fGeom->RelToAbsNumbering(iRelId2x2, iAbsId2x2);
	198	fGeom->RelToAbsNumbering(iRelIdnxn, iAbsIdnxn);
	199	fGeom->RelPosInAlice(iAbsId2x2, pos2x2);
	200	fGeom->RelPosInAlice(iAbsIdnxn, posnxn);
64df000d	201
	202	TArrayF triggerPosition(6);
	203	triggerPosition[0] = pos2x2(0) ;
	204	triggerPosition[1] = pos2x2(1) ;
	205	triggerPosition[2] = pos2x2(2) ;
	206	triggerPosition[3] = posnxn(0) ;
	207	triggerPosition[4] = posnxn(1) ;
	208	triggerPosition[5] = posnxn(2) ;
	209
	210	TArrayF triggerAmplitudes(4);
	211	triggerAmplitudes[0] = maxAmp2x2 ;
	212	triggerAmplitudes[1] = ampOutOfPatch2x2 ;
	213	triggerAmplitudes[2] = maxAmpnxn ;
	214	triggerAmplitudes[3] = ampOutOfPatchnxn ;
	215
	216	//esd->SetPHOSTriggerCells(triggerPosition);
	217	esd->AddPHOSTriggerPosition(triggerPosition);
	218	esd->AddPHOSTriggerAmplitudes(triggerAmplitudes);
	219
e68222ce	220
a5e00016	221	//########################################
	222	//############# Fill CaloCells ###########
	223	//########################################
	224
	225	Int_t nDigits = digitsArray->GetEntries();
	226	Int_t idignew = 0 ;
	227	AliDebug(1,Form("%d digits",nDigits));
	228
	229	const Int_t knEMC = fGeom->GetNModules()fGeom->GetNPhi()fGeom->GetNZ();
	230	AliESDCaloCells &phsCells = *(esd->GetPHOSCells());
	231	phsCells.CreateContainer(nDigits);
	232	phsCells.SetType(AliESDCaloCells::kPHOSCell);
	233
	234	// Add to CaloCells only EMC digits with non-zero energy
	235	for (Int_t idig = 0 ; idig < nDigits ; idig++) {
	236	const AliPHOSDigit * dig = (const AliPHOSDigit*)digitsArray->At(idig);
	237	if(dig->GetId() <= knEMC && dig->GetEnergy() > 0 ){
710634f4	238	//printf("i %d; id %d; amp %f; time %e\n",
710634f4	239	//idignew,dig->GetId(),dig->GetEnergy(), dig->GetTime());
a5e00016	240	phsCells.SetCell(idignew,dig->GetId(), dig->GetEnergy(), dig->GetTime());
	241	idignew++;
	242	}
e68222ce	243	}
a5e00016	244	phsCells.SetNumberOfCells(idignew);
a5e00016	245	phsCells.Sort();
e68222ce	246
a5e00016	247	//########################################
	248	//############## Fill CaloClusters #######
	249	//########################################
f4e0dd30	250
dfe0be07	251	for (Int_t recpart = 0 ; recpart < nOfRecParticles ; recpart++) {
a5e00016	252	AliPHOSRecParticle rp = dynamic_cast<AliPHOSRecParticle>(recParticles->At(recpart));
dfe0be07	253	if (Debug())
dfe0be07	254	rp->Print();
25ed816e	255	// Get track segment and EMC rec.point associated with this rec.particle
a5e00016	256	AliPHOSTrackSegment ts = static_cast<AliPHOSTrackSegment >(tsm->GetTrackSegments()
a5e00016	257	->At(rp->GetPHOSTSIndex()));
9a2cdbdf	258
a5e00016	259	AliPHOSEmcRecPoint emcRP = static_cast<AliPHOSEmcRecPoint >(emcRecPoints->At(ts->GetEmcIndex()));
25ed816e	260	AliESDCaloCluster *ec = new AliESDCaloCluster() ;
e68222ce	261
85c60a8e	262	Float_t xyz[3];
8013c27e	263	for (Int_t ixyz=0; ixyz<3; ixyz++)
8013c27e	264	xyz[ixyz] = rp->GetPos()[ixyz];
dd7ee508	265
dd7ee508	266	AliDebug(2,Form("Global position xyz=(%f,%f,%f)",xyz[0],xyz[1],xyz[2]));
77ea1c6f	267
a5e00016	268	// Create cell lists
	269
	270	Int_t cellMult = emcRP->GetDigitsMultiplicity();
	271	Int_t *digitsList = emcRP->GetDigitsList();
	272	Float_t *rpElist = emcRP->GetEnergiesList() ;
	273	UShort_t *absIdList = new UShort_t[cellMult];
	274	Double_t *fracList = new Double_t[cellMult];
	275
	276	for (Int_t iCell=0; iCell<cellMult; iCell++) {
	277	AliPHOSDigit digit = static_cast<AliPHOSDigit >(digitsArray->At(digitsList[iCell]));
	278	absIdList[iCell] = (UShort_t)(digit->GetId());
	279	if (digit->GetEnergy() > 0)
	280	fracList[iCell] = rpElist[iCell]/digit->GetEnergy();
	281	else
	282	fracList[iCell] = 0;
25ed816e	283	}
77ea1c6f	284
7592dfc4	285	//Primaries
7592dfc4	286	Int_t primMult = 0;
4dd59c4a	287	Int_t *primList = emcRP->GetPrimaries(primMult);
a5e00016	288
7592dfc4	289	// fills the ESDCaloCluster
8ada0ffe	290	ec->SetClusterType(AliESDCaloCluster::kPHOSCluster);
a5e00016	291	ec->SetPosition(xyz); //rec.point position in MARS
a5e00016	292	ec->SetE(rp->Energy()); //total particle energy
25ed816e	293	ec->SetClusterDisp(emcRP->GetDispersion()); //cluster dispersion
e68222ce	294	ec->SetPid(rp->GetPID()) ; //array of particle identification
25ed816e	295	ec->SetM02(emcRP->GetM2x()) ; //second moment M2x
	296	ec->SetM20(emcRP->GetM2z()) ; //second moment M2z
	297	ec->SetNExMax(emcRP->GetNExMax()); //number of local maxima
a5e00016	298	ec->SetEmcCpvDistance(ts->GetCpvDistance("r")); //Only radius, what about separate x,z????
25ed816e	299	ec->SetClusterChi2(-1); //not yet implemented
25ed816e	300	ec->SetM11(-1) ; //not yet implemented
3744a6cf	301
a5e00016	302	//Cells contributing to clusters
	303	ec->SetNCells(cellMult);
	304	ec->SetCellsAbsId(absIdList);
	305	ec->SetCellsAmplitudeFraction(fracList);
7592dfc4	306
3744a6cf	307	//Distance to the nearest bad crystal
3744a6cf	308	ec->SetDistanceToBadChannel(emcRP->GetDistanceToBadCrystal());
7592dfc4	309
7592dfc4	310	//Array of MC indeces
4dd59c4a	311	TArrayI arrayPrim(primMult,primList);
7592dfc4	312	ec->AddLabels(arrayPrim);
	313
	314	//Array of tracks uncomment when available in future
	315	//TArrayS arrayTrackMatched(1);// Only one track, temporal solution.
	316	//arrayTrackMatched[0]= (Short_t)(matchedTrack[iClust]);
	317	//ec->AddTracksMatched(arrayTrackMatched);
64df000d	318
dfe0be07	319	// add the track to the esd object
a5e00016	320
85c60a8e	321	esd->AddCaloCluster(ec);
9a2cdbdf	322	delete ec;
e68222ce	323	}
a5e00016	324	digitsArray ->Delete();
	325	delete digitsArray;
	326	emcRecPoints->Delete();
	327	delete emcRecPoints;
e68222ce	328	delete tsm;
e68222ce	329	delete pid;
dd7ee508	330	}
dd7ee508	331
e68222ce	332	//____________________________________________________________________________
d76c31f4	333	AliTracker* AliPHOSReconstructor::CreateTracker() const
dd7ee508	334	{
9a2cdbdf	335	// creates the PHOS tracker
	336	return new AliPHOSTracker();
	337	}
dd7ee508	338
1267d56d	339	//____________________________________________________________________________
9a2cdbdf	340	void AliPHOSReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digitsTree) const
	341	{
	342	// Converts raw data to
	343	// PHOS digits
	344	// Works on a single-event basis
dd7ee508	345
9a2cdbdf	346	rawReader->Reset() ;
dd7ee508	347
77ea1c6f	348	AliPHOSRawDecoder * dc ;
77ea1c6f	349
b3006690	350	const TObjArray* maps = AliPHOSRecoParamEmc::GetMappings();
	351	if(!maps) AliFatal("Cannot retrieve ALTRO mappings!!");
	352
	353	AliAltroMapping *mapping[4];
	354	for(Int_t i = 0; i < 4; i++) {
	355	mapping[i] = (AliAltroMapping*)maps->At(i);
	356	}
	357
77ea1c6f	358	if(strcmp(fgkRecoParamEmc->DecoderVersion(),"v1")==0)
b3006690	359	dc=new AliPHOSRawDecoderv1(rawReader,mapping);
710634f4	360	else
	361	if(strcmp(fgkRecoParamEmc->DecoderVersion(),"v2")==0)
	362	dc=new AliPHOSRawDecoderv2(rawReader,mapping);
	363	else
	364	dc=new AliPHOSRawDecoder(rawReader,mapping);
77ea1c6f	365
a7660309	366	dc->SetOldRCUFormat(fgkRecoParamEmc->IsOldRCUFormat());
86121c91	367
77ea1c6f	368	dc->SubtractPedestals(fgkRecoParamEmc->SubtractPedestals());
86121c91	369
9a2cdbdf	370	TClonesArray *digits = new TClonesArray("AliPHOSDigit",1);
	371	digits->SetName("DIGITS");
	372	Int_t bufsize = 32000;
	373	digitsTree->Branch("PHOS", &digits, bufsize);
dd7ee508	374
a28333c5	375	AliPHOSRawDigiProducer pr(fgkRecoParamEmc,fgkRecoParamCpv);
77ea1c6f	376	pr.MakeDigits(digits,dc);
	377
	378	delete dc ;
7592dfc4	379
9a2cdbdf	380	//!!!!for debug!!!
a28333c5	381	/*
9a2cdbdf	382	Int_t modMax=-111;
	383	Int_t colMax=-111;
	384	Int_t rowMax=-111;
	385	Float_t eMax=-333;
	386	//!!!for debug!!!
	387
	388	Int_t relId[4];
	389	for(Int_t iDigit=0; iDigit<digits->GetEntries(); iDigit++) {
	390	AliPHOSDigit* digit = (AliPHOSDigit*)digits->At(iDigit);
	391	if(digit->GetEnergy()>eMax) {
	392	fGeom->AbsToRelNumbering(digit->GetId(),relId);
	393	eMax=digit->GetEnergy();
	394	modMax=relId[0];
	395	rowMax=relId[2];
	396	colMax=relId[3];
dd7ee508	397	}
	398	}
	399
9a2cdbdf	400	AliDebug(1,Form("Digit with max. energy: modMax %d colMax %d rowMax %d eMax %f\n\n",
9a2cdbdf	401	modMax,colMax,rowMax,eMax));
a28333c5	402	*/
9a2cdbdf	403	digitsTree->Fill();
e68222ce	404	digits->Delete();
e68222ce	405	delete digits;
d7d3b67b	406	}