splits clusters on the basis of exponential distance from the bump centroid More...

#include <PndEmcExpClusterSplitter.h>

Inheritance diagram for PndEmcExpClusterSplitter:

Public Member Functions
	PndEmcExpClusterSplitter (Int_t verbose=0)

virtual	~PndEmcExpClusterSplitter ()

virtual InitStatus	Init ()
	Init Task. More...

virtual void	Exec (Option_t *opt)
	Runs the task. More...

virtual void	FinishTask ()
	Called at end of task. More...

void	SetStorageOfData (Bool_t p=kTRUE)

PndEmcBump *	AddBump ()
	Adds a new PndEmcBump to fBumpArray and returns it. More...

PndEmcSharedDigi *	AddSharedDigi (PndEmcDigi *, Double_t weight)
	Adds a new PndEmcSharedDigi to fSharedDigiArray and returns it. More...

void	SetPersistency (Bool_t val=kTRUE)

Bool_t	GetPersistency ()

Protected Member Functions
virtual void	SetParContainers ()

Private Member Functions
	PndEmcExpClusterSplitter (const PndEmcExpClusterSplitter &)

PndEmcExpClusterSplitter &	operator= (const PndEmcExpClusterSplitter &)

	ClassDef (PndEmcExpClusterSplitter, 2)

Private Attributes
TClonesArray *	fDigiArray

TClonesArray *	fClusterArray

TClonesArray *	fBumpArray

TClonesArray *	fSharedDigiArray

PndEmcGeoPar *	fGeoPar
	Geometry parameter container. More...

PndEmcDigiPar *	fDigiPar
	Digitisation parameter container. More...

PndEmcRecoPar *	fRecoPar
	Reconstruction parameter container. More...

std::vector< Double_t >	fClusterPosParam

Double_t	fMoliereRadius

Double_t	fMoliereRadiusShashlyk

Double_t	fExponentialConstant

Int_t	fMaxIterations

Double_t	fCentroidShift

Int_t	fMaxBumps

Double_t	fMinDigiEnergy

PndEmcDigiCalibrator	digiCalibrator

Int_t	HowManyDidis

Detailed Description

splits clusters on the basis of exponential distance from the bump centroid

Definition at line 57 of file PndEmcExpClusterSplitter.h.

Constructor & Destructor Documentation

PndEmcExpClusterSplitter::PndEmcExpClusterSplitter ( Int_t verbose = 0 )

Definition at line 71 of file PndEmcExpClusterSplitter.cxx.

References fClusterPosParam, and PndPersistencyTask::SetPersistency().

                                                                :PndPersistencyTask("PndEmcExpClusterSplitter", verbose),
 fDigiArray(0), fClusterArray(0), fBumpArray(0), fSharedDigiArray(0), fGeoPar(new PndEmcGeoPar()), fDigiPar(new PndEmcDigiPar()), fRecoPar(new PndEmcRecoPar()), fClusterPosParam(), fMoliereRadius(0), fMoliereRadiusShashlyk(0), fExponentialConstant(0), fMaxIterations(0), fCentroidShift(0), fMaxBumps(0), fMinDigiEnergy(0)
 {
   fClusterPosParam.clear();
   SetPersistency(kTRUE);
 }

PndEmcExpClusterSplitter::~PndEmcExpClusterSplitter ( )

virtual

Definition at line 82 of file PndEmcExpClusterSplitter.cxx.

 {
 //      delete fGeoPar;
 //      delete fDigiPar;
 //      delete fRecoPar;
 }

PndEmcExpClusterSplitter::PndEmcExpClusterSplitter ( const PndEmcExpClusterSplitter & )

private

Member Function Documentation

PndEmcBump * PndEmcExpClusterSplitter::AddBump ( )

Adds a new PndEmcBump to fBumpArray and returns it.

Returns: PndEmcBump*

Definition at line 566 of file PndEmcExpClusterSplitter.cxx.

References fBumpArray.

Referenced by Exec().

 {
         TClonesArray& clref = *fBumpArray;
         Int_t size = clref.GetEntriesFast();
         return new(clref[size]) PndEmcBump();
 }

PndEmcSharedDigi * PndEmcExpClusterSplitter::AddSharedDigi	(	PndEmcDigi *	digi,
		Double_t	weight
	)

Adds a new PndEmcSharedDigi to fSharedDigiArray and returns it.

Parameters

digi	Digi which is shared
weight	Weight of digi in this shared digi

Returns: PndEmcSharedDigi*

Definition at line 580 of file PndEmcExpClusterSplitter.cxx.

References fSharedDigiArray.

Referenced by Exec().

                                                                                           {
         TClonesArray& clref = *fSharedDigiArray;
         Int_t size = clref.GetEntriesFast();
         return new(clref[size]) PndEmcSharedDigi(*digi, weight);
 }

PndEmcExpClusterSplitter::ClassDef	(	PndEmcExpClusterSplitter	,
		2
	)

private

void PndEmcExpClusterSplitter::Exec ( Option_t * opt )

virtual

Runs the task.

The algorithm is as follows: We will index each bump by its maximum digi's PndEmcTwoCoordIndex. We will set up a list of bump centroids which to start with will be synonymous with the location of the maxima. We then apportion a weight to each digi, according to its distance from the centroids. We then construct the bumps according to these weights, which will presumably give a different set of centroids. This is repeated until the centroids are static within tolerance, or we reach the maximum number of iterations.

Parameters

opt unused

Returns: void

Definition at line 189 of file PndEmcExpClusterSplitter.cxx.

References PndEmcXClMoments::AbsZernikeMoment(), AddBump(), PndEmcCluster::addDigi(), AddSharedDigi(), PndEmcDigiCalibrator::CalibrationEvtTimeByDigi(), digiCalibrator, PndEmcCluster::DigiList(), Double_t, PndEmcClusterProperties::Energy(), exp(), fBumpArray, fCentroidShift, fClusterArray, fClusterPosParam, fDigiArray, PndEmcDigi::fEvtNo, fExponentialConstant, fMaxBumps, fMinDigiEnergy, fMoliereRadius, fMoliereRadiusShashlyk, fRecoPar, fSharedDigiArray, fVerbose, PndEmcRecoPar::GetEmcClusterPosMethod(), PndEmcSharedDigi::GetEnergy(), PndEmcDigi::GetEnergy(), PndEmcDigi::GetModule(), PndEmcMapper::GetTCI(), PndEmcDigi::GetTCI(), PndEmcDigiCalibrator::GetTimeResolutionOfDigi(), PndEmcCluster::GetTrackEntering(), HowManyDidis, i, PndEmcTwoCoordIndex::Index(), PndEmcMapper::Instance(), PndEmcCluster::IsEnergyValid(), PndEmcCluster::IsPositionValid(), PndEmcXClMoments::Lat(), PndEmcCluster::LocalMaxMap(), PndEmcBump::MadeFrom(), map, PndEmcCluster::MemberDigiMap(), PndEmcCluster::NumberOfDigis(), PndEmcCluster::SetEnergy(), PndEmcBump::SetEventNo(), PndEmcCluster::SetLatMom(), PndEmcCluster::SetNBumps(), PndEmcCluster::SetPosition(), PndEmcCluster::SetZ20(), PndEmcCluster::SetZ53(), PndEmcClusterProperties::Where(), PndEmcDigi::where(), PndEmcTwoCoordIndex::XCoord(), and PndEmcTwoCoordIndex::YCoord().

 {
 
         PndEmcMapper *fEmcMap=PndEmcMapper::Instance();
         // Reset output array
         if ( ! fBumpArray ) Fatal("Exec", "No Bump Array");
         fBumpArray->Delete();
         if ( ! fSharedDigiArray ) Fatal("Exec", "No Shared Digi Array");
         fSharedDigiArray->Delete();
         
         int nClusters = fClusterArray->GetEntriesFast();
         
         for (Int_t iCluster = 0; iCluster < nClusters; iCluster++){
                 PndEmcCluster* theCluster = (PndEmcCluster*) fClusterArray->At(iCluster);
 
                 int numberOfBumps = -1;
                 numberOfBumps = (theCluster->LocalMaxMap()).size();
         
                 if (numberOfBumps<=1 || numberOfBumps >= fMaxBumps){
                         // Limit the max number of bumps in the cluster to 8 (default)
                         // in this case, we clearly have a cluster, but no bumps to speak of.  
                         // Make 1 bump with weights all equal to 1
                 
                         std::map<Int_t, Int_t>::const_iterator theDigiIterator;
                         
                         PndEmcBump* theNewBump = AddBump(); 
                         theNewBump->MadeFrom(iCluster);
                         theNewBump->SetLink(FairLink("EmcCluster", iCluster));
 
                         for(theDigiIterator = theCluster->MemberDigiMap().begin();
                         theDigiIterator != theCluster->MemberDigiMap().end(); ++theDigiIterator){
                                 PndEmcDigi *theDigi = (PndEmcDigi *) fDigiArray->At(theDigiIterator->second);
                                 AddSharedDigi(theDigi, 1.0); // PndEmcSharedDigi* sharedDigi= //[R.K.03/2017] unused variable
                                 Int_t iSharedDigi=fSharedDigiArray->GetEntriesFast()-1;
                                 theNewBump->addDigi(fSharedDigiArray,iSharedDigi);
                         }
                 }
                 else {
                         std::map<Int_t,Int_t> theMaximaDigis=theCluster->LocalMaxMap();
                         std::map<Int_t,Int_t>::iterator theMaximaDigisIterator;
                         std::map<PndEmcTwoCoordIndex*, TVector3*> theCentroidPoints;
                         std::map<PndEmcTwoCoordIndex*, TVector3*> theMaximaPoints;
                         std::map<PndEmcTwoCoordIndex*, PndEmcBump*> theIndexedBumps;
                         std::map<PndEmcTwoCoordIndex*, TVector3*> theAllDigiPoints;
                         
                         std::map<Int_t, Int_t> theDigiDict = theCluster->MemberDigiMap();
                         
                         double totalEnergy=0;
                         
                         for(theMaximaDigisIterator = theMaximaDigis.begin(); 
                         theMaximaDigisIterator != theMaximaDigis.end();
                         ++theMaximaDigisIterator){
                                 PndEmcDigi *theMaxDigi = (PndEmcDigi *) fDigiArray->At(theMaximaDigisIterator->second);
                                 
                                 Int_t detId = theMaximaDigisIterator->first;
                                 PndEmcTwoCoordIndex *theTCI = fEmcMap->GetTCI(detId);
                                 
                                 totalEnergy += theMaxDigi->GetEnergy();
                                 
                                 TVector3 *digiLocation = new TVector3(theMaxDigi->where());
                                 TVector3 *sameLocation = new TVector3(theMaxDigi->where());
                                 
                                 theMaximaPoints.insert(std::map<PndEmcTwoCoordIndex*, TVector3*>::value_type(theTCI, digiLocation));
                                 theCentroidPoints.insert(std::map<PndEmcTwoCoordIndex*, TVector3*>::value_type(theTCI, sameLocation));
                         }
         
                         std::map<Int_t,Int_t>::iterator theAllDigisIterator;            
                         // This loop works out the location of all the digis in the cluster 
                 
                         for( theAllDigisIterator = theDigiDict.begin(); theAllDigisIterator != theDigiDict.end();
                                         ++theAllDigisIterator){
                                 Int_t detId = theAllDigisIterator->first;
                                 PndEmcTwoCoordIndex *theTCI = fEmcMap->GetTCI(detId);
                                 PndEmcDigi *theDigi = (PndEmcDigi *) fDigiArray->At(theAllDigisIterator->second);
                                 TVector3 *digiLocation = new TVector3(theDigi->where());
                                 theAllDigiPoints.insert(std::map<PndEmcTwoCoordIndex*, TVector3*>::value_type( theTCI, digiLocation ));
                         }
 
                         theMaximaDigisIterator = theMaximaDigis.begin();
 
                         // Now we can create the EmcBumps
 
                         // The algorithm is as follows: We will index each bump by its
                         // maximum digi's PndEmcTwoCoordIndex.  We will set up a list of
                         // bump centroids which to start with will be synonymous with the
                         // location of the maxima.  We then apportion a weight to each
                         // digi, according to its distance from the centroids.  We then
                         // construct the bumps according to these weights, which will
                         // presumably give a different set of centroids.  This is repeated
                         // until the centroids are static within tolerance, or we reach
                         // the maximum number of iterations.
 
                         Int_t iterations = 0;
 
                         Double_t averageCentroidShift;
 
                         do {
                                 if (fVerbose>=3){
                                         std::cout<<"iteration No "<<iterations<<std::endl;
                                 }
                                 averageCentroidShift=0.0;
 
                                 // First clean up the old bumps
                                 std::map<PndEmcTwoCoordIndex*, PndEmcBump*>::iterator theBumpKiller = theIndexedBumps.begin();
                                 while(theBumpKiller != theIndexedBumps.end()){
                                         PndEmcBump *theBump = theBumpKiller->second;
                                         delete theBump;
                                         ++theBumpKiller;
                                 }
                                 theIndexedBumps.clear();
 
                                 // Then loop over all the maxima and assign weights accordingly
                                 for (theMaximaDigisIterator = theMaximaDigis.begin();
                                                 theMaximaDigisIterator != theMaximaDigis.end();
                                                 ++theMaximaDigisIterator) {
                                         Int_t detId = theMaximaDigisIterator->first;
                                         PndEmcTwoCoordIndex *theCurrentMaximaTCI = fEmcMap->GetTCI(detId);
 
                                         if (fVerbose>=3){
                                                 std::cout<<"***************** current maximum: theta = "<<theCurrentMaximaTCI->XCoord()
                                                         <<", phi = "<<theCurrentMaximaTCI->YCoord()<<"*********"<<std::endl;
                                         }
 
                                         // Create the bump which will correspond to this digi maxima
                                         PndEmcBump* theNewBump = new PndEmcBump();
                                         theNewBump->MadeFrom(iCluster);
                                         theIndexedBumps.insert(std::map<PndEmcTwoCoordIndex*,
                                                         PndEmcBump*>::value_type(theCurrentMaximaTCI, theNewBump));
 
 
                                         // Now we will look over all the digis and add each of them
                                         // to this Bump with an appropriate weight
 
                                         for (theAllDigisIterator = theDigiDict.begin();
                                                         theAllDigisIterator != theDigiDict.end();++theAllDigisIterator) {
                                                 PndEmcDigi *theCurrentDigi = (PndEmcDigi *) fDigiArray->At(theAllDigisIterator->second);
                                                 PndEmcTwoCoordIndex *theCurrentTCI = theCurrentDigi->GetTCI();
 
                                                 Double_t weight;
 
                                                 // We are on the first pass and the digi is not a local max, or we are not on the
                                                 // first pass.   Assign a weight according to the distance from the centroid position.
 
                                                 Double_t myEnergy = 0;
                                                 Double_t myDistance = 0;
 
                                                 // Now share the digi out according to its distance from the maxima, 
                                                 // and the maxima energies
 
                                                 Double_t totalDistanceEnergy=0;
 
                                                 //Moliere Radius for Shashlyk is different
                                                 Double_t MoliereRadius;
                                                 if(theCurrentDigi->GetModule() == 5)
                                                         MoliereRadius = fMoliereRadiusShashlyk;
                                                 else
                                                         MoliereRadius = fMoliereRadius;
 
                                                 std::map<PndEmcTwoCoordIndex*,TVector3*>::iterator theMaxPointsIterator;
 
                                                 for(theMaxPointsIterator = theCentroidPoints.begin(); theMaxPointsIterator != theCentroidPoints.end();++theMaxPointsIterator){
                                                         PndEmcTwoCoordIndex *theMaxPointsTCI =theMaxPointsIterator->first; 
 
                                                         TVector3 *theMaxPoint = theMaxPointsIterator->second;
                                                         TVector3 *theCurrentDigiPoint = theAllDigiPoints.find(theCurrentTCI)->second;
 
                                                         Double_t theDistance;
 
                                                         // This next bit just checks to see if the maxima point in
                                                         // hand is the same as the crystal from which we are
                                                         // trying to find distance - just an FP trap really.
 
                                                         if ((*theCurrentTCI)==(*theMaxPointsTCI)){
                                                                 theDistance=0.0;
                                                         } else {
                                                                 TVector3 distance( *theMaxPoint - *theCurrentDigiPoint);
 
                                                                 theDistance = distance.Mag();
                                                         }
 
                                                         if (*theCurrentMaximaTCI == *(theMaxPointsTCI)){
                                                                 // i.e. the maximum we are trying to find the distance from is 
                                                                 // the one for which we are currently trying to make a bump
                                                                 myDistance = theDistance;
                                                                 Int_t iCurentMaxDigi = (theDigiDict.find(theMaxPointsTCI->Index()))->second;
                                                                 myEnergy = ((PndEmcDigi *) fDigiArray->At(iCurentMaxDigi))->GetEnergy();
                                                         }
 
                                                         Int_t iMaxPoint = (theDigiDict.find(theMaxPointsTCI->Index()))->second;
                                                         totalDistanceEnergy += ((PndEmcDigi *) fDigiArray->At(iMaxPoint))->GetEnergy() *
                                                                 exp(-fExponentialConstant * theDistance/MoliereRadius);
                                                 }
 
                                                 if(totalDistanceEnergy > 0.0)
                                                         weight = myEnergy*exp(-fExponentialConstant* 
                                                                         myDistance/MoliereRadius) / ( totalDistanceEnergy);
                                                 else 
                                                         weight=0;
 
                                                 if (fVerbose>=3){
                                                         std::cout<<"\t digi theta = "<<theCurrentDigi->GetTCI()->XCoord()
                                                                 <<", phi = "<<theCurrentDigi->GetTCI()->YCoord()<<std::endl;
                                                         std::cout<<"energy = "<<theCurrentDigi->GetEnergy()<<", weight = "<< weight<<std::endl;
                                                 }
                                                 PndEmcSharedDigi* sharedDigi= AddSharedDigi( theCurrentDigi, weight );
 
                                                 if (fVerbose>=3){
                                                         std::cout<<"shared digi energy = "<<sharedDigi->GetEnergy()<<std::endl;
                                                 }
 
                                                 Int_t iSharedDigi=fSharedDigiArray->GetEntriesFast()-1;
                                                 if( sharedDigi->GetEnergy() > fMinDigiEnergy){
                                                         theNewBump->addDigi( fSharedDigiArray, iSharedDigi );
                                                 } else {
                                                         fSharedDigiArray->RemoveAt(iSharedDigi);
                                                         fSharedDigiArray->Compress();
                                                 }
                                         }
 
                                         // Compute the shift of the centroid we have just calculated
                                         TVector3 *theOldCentroid = theCentroidPoints.find(theCurrentMaximaTCI)->second;
 
                                         PndEmcClusterProperties clusterProperties(*theNewBump, fSharedDigiArray);
 
                                         TVector3 newbumppos =  clusterProperties.Where(fRecoPar->GetEmcClusterPosMethod(), fClusterPosParam);
                                         theNewBump->SetPosition(newbumppos);
                                         TVector3 centroidShift(*theOldCentroid - newbumppos);
                                         averageCentroidShift+=centroidShift.Mag();
                                 }
 
                                 averageCentroidShift/=(Double_t)numberOfBumps;
 
                                 // Put the new centroids in the list of centroid points,
                                 // remembering to delete the old ones.
                                 std::map<PndEmcTwoCoordIndex*, TVector3*>::iterator 
                                         theCentroidPointsIterator = theCentroidPoints.begin();
                                 for(theCentroidPointsIterator = theCentroidPoints.begin();
                                                 theCentroidPointsIterator != theCentroidPoints.end(); 
                                                 ++theCentroidPointsIterator) {
                                         delete theCentroidPointsIterator->second;
                                 }
                                 theCentroidPoints.clear();
 
                                 std::map<PndEmcTwoCoordIndex*, PndEmcBump*>::iterator theIndexedBumpsIterator;
                                 for(theIndexedBumpsIterator = theIndexedBumps.begin();
                                                 theIndexedBumpsIterator != theIndexedBumps.end(); ++theIndexedBumpsIterator){
                                         TVector3 *theNewCentroid = new TVector3((theIndexedBumpsIterator->second)->where());
                                         theCentroidPoints.insert(std::map<PndEmcTwoCoordIndex*, TVector3*>::value_type(theIndexedBumpsIterator->first,theNewCentroid));
                                 }
 
                                 iterations++;
 
                         } while (iterations < fMaxIterations && averageCentroidShift > fCentroidShift);
                         //End of do loop
 
                         // Finally append the new bumps to the TClonesArray.
                         std::map<PndEmcTwoCoordIndex*, PndEmcBump*>::iterator theBumpsIterator;
                         PndEmcBump *theBump;
 
                         for(theBumpsIterator = theIndexedBumps.begin(); theBumpsIterator != theIndexedBumps.end();++theBumpsIterator){
                                 theBump = theBumpsIterator->second;
                                 Int_t size_ba = fBumpArray->GetEntriesFast();
                                 PndEmcBump* theNextBump = new((*fBumpArray)[size_ba]) PndEmcBump(*(theBump));
                                 if (fVerbose>0)
                                         std::cout << "Bump Created!" << std::endl;
                                 theNextBump->SetInsertHistory(kFALSE);
                                 theNextBump->SetLink(FairLink("EmcCluster", iCluster));
                                 PndEmcCluster* myCluster = (PndEmcCluster*)fClusterArray->At(iCluster);
                                 theNextBump->AddLinks(myCluster->GetTrackEntering());
                                 theNextBump->SetTimeStamp(myCluster->GetTimeStamp());
                                 theNextBump->SetTimeStampError(myCluster->GetTimeStampError());
                         }
 
                         std::map<PndEmcTwoCoordIndex*,TVector3*>::iterator theGrimReaper = theMaximaPoints.begin();
                         for(theGrimReaper = theMaximaPoints.begin();
                                         theGrimReaper != theMaximaPoints.end();
                                         ++theGrimReaper){
                                 delete theGrimReaper->second;
                         }
                         theMaximaPoints.clear();
 
                         for(theGrimReaper  = theAllDigiPoints.begin();
                                         theGrimReaper != theAllDigiPoints.end();
                                         ++theGrimReaper){
                                 delete theGrimReaper->second;
                         }
                         theAllDigiPoints.clear();
 
                         for(theGrimReaper = theCentroidPoints.begin(); 
                                         theGrimReaper != theCentroidPoints.end();
                                         ++theGrimReaper){
                                 delete theGrimReaper->second;
                         }
                         theCentroidPoints.clear();
                 }
 
                 Int_t nBumps = (theCluster->LocalMaxMap()).size();
                 theCluster->SetNBumps(nBumps);
         }
 
         // At that moment internal state fEnergy and fWhere of Clusters are
         // not initialized, the following make it possible to see energy and
         // position from output root file
         Int_t nBump = fBumpArray->GetEntriesFast();
 
         Double_t  fTimeError; //CalibTimeOfaDigi, //[R.K.03/2017] unused variable
         Double_t WeightedFactor1(0.), NormWeightedFactor1(0.), AverageTime1(0.);
         //Double_t WeightedFactor2(0.), NormWeightedFactor2(0.), AverageTime2(0.);
         //Double_t WeightedFactor3(0.), NormWeightedFactor3(0.), AverageTime3(0.);
         for (Int_t i=0; i<nBump; i++){
                 PndEmcBump *tmpbump = (PndEmcBump*) fBumpArray->At(i);
                 PndEmcClusterProperties clusterProperties(*tmpbump, fSharedDigiArray);
 
                 if (!tmpbump->IsEnergyValid())
                         tmpbump->SetEnergy(clusterProperties.Energy());
                 if (!tmpbump->IsPositionValid())
                         tmpbump->SetPosition(clusterProperties.Where(fRecoPar->GetEmcClusterPosMethod(), fClusterPosParam));
                 PndEmcXClMoments xClMoments(*tmpbump, fSharedDigiArray);
                 tmpbump->SetZ20(xClMoments.AbsZernikeMoment(2, 0, 15));
                 tmpbump->SetZ53(xClMoments.AbsZernikeMoment(5, 3, 15));
                 tmpbump->SetLatMom(xClMoments.Lat());
 
                 //for time information
                 WeightedFactor1 = 0.;//= WeightedFactor2 = WeightedFactor3 = 0.;
                 AverageTime1 = 0.;//AverageTime2 = AverageTime3 = 0.;
                 Double_t fMaxDigiEnergy = -1.;
                 const std::vector<Int_t>& listOfDigi = tmpbump->DigiList();
                 for(size_t id=0;id <listOfDigi.size();++id){
                         PndEmcDigi* theDigi = (PndEmcDigi*)fSharedDigiArray->At(listOfDigi[id]);
                         //CalibTimeOfaDigi = digiCalibrator.CalibrationEvtTimeByDigi(theDigi, kFALSE);
                         fTimeError = digiCalibrator.GetTimeResolutionOfDigi(theDigi);
                         WeightedFactor1 += 1./fTimeError/fTimeError;
                         //WeightedFactor2 += theDigi->GetEnergy()*1./fTimeError/fTimeError;
                         //WeightedFactor3 += theDigi->GetEnergy();
                         if(theDigi->GetEnergy() > fMaxDigiEnergy){
                                 fMaxDigiEnergy = theDigi->GetEnergy();
                                 tmpbump->SetEventNo(theDigi->fEvtNo);
                                 //tmpbump->SetTimeStamp1(theDigi->GetTimeStamp());
                                 //tmpbump->fSeedPosition = theDigi->where();
                         }
                 }
                 for(size_t id=0;id <listOfDigi.size();++id){
                         PndEmcDigi* theDigi = (PndEmcDigi*)fSharedDigiArray->At(listOfDigi[id]);
                         digiCalibrator.CalibrationEvtTimeByDigi(theDigi, kFALSE); //CalibTimeOfaDigi =  //[R.K.03/2017] unused variable
                         fTimeError = digiCalibrator.GetTimeResolutionOfDigi(theDigi);
                         NormWeightedFactor1 = 1./fTimeError/fTimeError;
                         NormWeightedFactor1 /= WeightedFactor1;
                         AverageTime1 += NormWeightedFactor1*theDigi->GetTimeStamp();
 
                         //NormWeightedFactor2 = theDigi->GetEnergy()*1./fTimeError/fTimeError;
                         //NormWeightedFactor2 /= WeightedFactor2;
                         //AverageTime2 += NormWeightedFactor2*theDigi->GetTimeStamp();
 
                         //NormWeightedFactor3 = theDigi->GetEnergy();
                         //NormWeightedFactor3 /= WeightedFactor3;
                         //AverageTime3 += NormWeightedFactor3*theDigi->GetTimeStamp();
                 }
                 tmpbump->SetTimeStamp(AverageTime1);
                 HowManyDidis += tmpbump->NumberOfDigis();
                 //tmpbump->SetTimeStamp1(AverageTime1);
                 //tmpbump->SetTimeStamp2(AverageTime2);
                 //tmpbump->SetTimeStamp3(AverageTime3);
                 //end for time information
         }
 
         if (fVerbose>=1){
                 std::cout<<"PndEmcExpClusterSplitter:: Number of clusters = "<<nClusters<<std::endl;
                 std::cout<<"PndEmcExpClusterSplitter:: Number of bumps = "<<nBump<<std::endl;
         }
 
 }

void PndEmcExpClusterSplitter::FinishTask ( )

virtual

Called at end of task.

Outputs the number of digis read.

Returns: void

Definition at line 593 of file PndEmcExpClusterSplitter.cxx.

References HowManyDidis.

 {
         cout<<"================================================="<<endl;
         cout<<"PndEmcExpClusterSplitter::FinishTask"<<endl;
         cout<<"================================================="<<endl;
         cout<<"read digis #"<<HowManyDidis<<endl;
 }

Bool_t PndPersistencyTask::GetPersistency ( )

inlineinherited

Definition at line 32 of file PndPersistencyTask.h.

References PndPersistencyTask::fPersistency.

32 { return fPersistency; }

PndPersistencyTask::fPersistency

Bool_t fPersistency

Definition: PndPersistencyTask.h:35

InitStatus PndEmcExpClusterSplitter::Init ( )

virtual

Init Task.

Returns: InitStatus

Return values

kSUCCESS success

Definition at line 95 of file PndEmcExpClusterSplitter.cxx.

                                           {
   
   // Get RootManager
   FairRootManager* ioman = FairRootManager::Instance();
   if ( ! ioman ){
     cout << "-E- PndEmcExpClusterSplitter::Init: "
          << "RootManager not instantiated!" << endl;
     return kFATAL;
   }
         
   // Geometry loading
   fGeoPar->InitEmcMapper();
   PndEmcStructure::Instance();
         
         // Get input array
         if(FairRunAna::Instance()->IsTimeStamp()) {
                 fDigiArray = (TClonesArray*) ioman->GetObject("EmcDigiClusterBase");
         } else {
                 fDigiArray = (TClonesArray*) ioman->GetObject("EmcDigi");
         }
         if ( ! fDigiArray ) {
                 cout << "-W- PndEmcExpClusterSplitter::Init: "
                 << "No PndEmcDigi array!" << endl;
                 return kERROR;
         }
   
   fClusterArray = dynamic_cast<TClonesArray *> (ioman->GetObject("EmcCluster"));
   if ( ! fClusterArray ) {
     cout << "-W- PndEmcExpClusterSplitter::Init: "
          << "No PndEmcCluster array!" << endl;
     return kERROR;
   }
   
   fMoliereRadius=fRecoPar->GetMoliereRadius();// Mr in cm
   fMoliereRadiusShashlyk=fRecoPar->GetMoliereRadiusShashlyk();// Mr in cm
   fExponentialConstant=fRecoPar->GetExponentialConstant();
   // Energy fall off with distance from centre of cluster is
   // exp(-a*dist/Mr) Mr is the moliere radius.  dist is distance from
   // centre in cm, a is the above parameter.  The optimized value for
   // logarithimc cluster positioning is 2.5 For linear cluster
   // positioning a should be set to 1.5
   
   fMaxIterations=fRecoPar->GetMaxIterations();
   // Set the max number of iterations that the splitting algorithm is allowed 
   // to do before it decides that enough is enough.
   
   fCentroidShift=fRecoPar->GetCentroidShift();
   // Set the tolerance level to which it is said that a bump is 
   // said not to have moved since the last iteration.  
   // The default is a millimetre.
   fMaxBumps=fRecoPar->GetMaxBumps();
   // Set an upper limit on the number of bumps allowed in a cluster if it 
   // is to be split.
   
   fMinDigiEnergy=fRecoPar->GetMinDigiEnergy();
   // Set minimum SharedDigi energy to 20keV.
         
         if (!strcmp(fRecoPar->GetEmcClusterPosMethod(),"lilo"))
         {
                 cout<<"Lilo cluster position method"<<endl;
                 fClusterPosParam.push_back(fRecoPar->GetOffsetParmA());
                 fClusterPosParam.push_back(fRecoPar->GetOffsetParmB());
                 fClusterPosParam.push_back(fRecoPar->GetOffsetParmC());
         }       
 
   // Create and register output array
   fBumpArray = new TClonesArray("PndEmcBump");
   ioman->Register("EmcBump","Emc",fBumpArray,GetPersistency());
   
   fSharedDigiArray = new TClonesArray("PndEmcSharedDigi");
   ioman->Register("EmcSharedDigi","Emc",fSharedDigiArray,GetPersistency());
 
         HowManyDidis = 0;
 
   cout << "-I- PndEmcExpClusterSplitter: Intialization successfull" << endl;
   return kSUCCESS;
 }

PndEmcExpClusterSplitter& PndEmcExpClusterSplitter::operator= ( const PndEmcExpClusterSplitter & )

private

void PndEmcExpClusterSplitter::SetParContainers ( )

protectedvirtual

Get parameter containers

Definition at line 601 of file PndEmcExpClusterSplitter.cxx.

References fDigiPar, fGeoPar, fRecoPar, and run.

                                                 {
 
         // Get run and runtime database
         FairRun* run = FairRun::Instance();
         if ( ! run ) Fatal("SetParContainers", "No analysis run");
 
         FairRuntimeDb* db = run->GetRuntimeDb();
         if ( ! db ) Fatal("SetParContainers", "No runtime database");
         // Get Emc digitisation parameter container
         fGeoPar = (PndEmcGeoPar*) db->getContainer("PndEmcGeoPar");
         // Get Emc digitisation parameter container
         fDigiPar = (PndEmcDigiPar*) db->getContainer("PndEmcDigiPar");
         // Get Emc reconstruction parameter container
         fRecoPar = (PndEmcRecoPar*) db->getContainer("PndEmcRecoPar");
 }