// --- Standard library ---
// --- AliRoot header files ---
-#include "AliPHOSCalibrationDB.h"
#include "AliPHOSClusterizerv1.h"
#include "AliPHOSCpvRecPoint.h"
#include "AliPHOSDigit.h"
#include "AliPHOS.h"
#include "AliPHOSGetter.h"
#include "AliRun.h"
+#include "AliPHOSCalibrManager.h"
+#include "AliPHOSCalibrationData.h"
ClassImp(AliPHOSClusterizerv1)
//____________________________________________________________________________
AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
{
- // dtor
+ delete fPedestals ;
+ delete fGains ;
+
fSplitFile = 0 ;
}
//____________________________________________________________________________
Float_t AliPHOSClusterizerv1::Calibrate(Int_t amp, Int_t absId) const
-{ //To be replased later by the method, reading individual parameters from the database
- if(fCalibrationDB)
- return fCalibrationDB->Calibrate(amp,absId) ;
+{
+ if(fPedestals || fGains ){ //use calibration data
+ if(!fPedestals || !fGains ){
+ Error("Calibrate","Either Pedestals of Gains not set!") ;
+ return 0 ;
+ }
+ Float_t en=(amp - fPedestals->Data(absId))*fGains->Data(absId) ;
+ if(en>0)
+ return en ;
+ else
+ return 0. ;
+ }
else{ //simulation
if(absId <= fEmcCrystals) //calibrate as EMC
return fADCpedestalEmc + amp*fADCchanelEmc ;
void AliPHOSClusterizerv1::Exec(Option_t * option)
{
// Steering method
-
if( strcmp(GetName(), "")== 0 )
Init() ;
AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
if(gime->BranchExists("RecPoints"))
return ;
+
+ //test, if this is real data:
+ TBranch * eh = gAlice->TreeE()->GetBranch("AliPHOSBeamTestEvent") ;
+ if(eh){ //Read CalibrationData
+ AliPHOSCalibrManager * cmngr=AliPHOSCalibrManager::GetInstance() ;
+ if(!cmngr){ //not defined yet
+ cmngr=AliPHOSCalibrManager::GetInstance("PHOSBTCalibration.root") ;
+ }
+ if(fPedestals)
+ delete fPedestals ;
+ fPedestals = new AliPHOSCalibrationData("Pedestals",fCalibrVersion) ;
+ cmngr->ReadFromRoot(*fPedestals,fCalibrRun) ;
+ if(fGains)
+ delete fGains ;
+ fGains = new AliPHOSCalibrationData("Gains",fCalibrVersion) ;
+ cmngr->ReadFromRoot(*fGains,fCalibrRun) ;
+ }
+
Int_t nevents = gime->MaxEvent() ;
Int_t ievent ;
for(ievent = 0; ievent < nevents; ievent++){
-
gime->Event(ievent,"D") ;
-
+ Int_t pattern = gime->EventPattern() ;
+ if(pattern){
+ printf("pattern %d \n",pattern) ;
+ if(!fPatterns){
+ Error("Exec","Patterns for reconstruction of events are not set, use SetPatterns() ") ;
+ return ;
+ }
+ Bool_t skip = kTRUE ;
+ Int_t npat = fPatterns->GetSize() ;
+ for(Int_t i = 0; i<npat;i++)
+ if(pattern==fPatterns->At(i)){
+ skip = kFALSE ;
+ break ;
+ }
+ if(skip)
+ continue ;
+ }
+
if(ievent == 0)
GetCalibrationParameters() ;
fADCpedestalCpv = dig->GetCPVpedestal() ;
}
else{
- fCalibrationDB = gime->CalibrationDB();
+ ;
+// fCalibrationDB = gime->CalibrationDB();
}
}
fEmcTimeGate = 1.e-8 ;
fToUnfold = kTRUE ;
-
+
+ fPurifyThreshold = 0.012 ;
+
+ fPedestals = 0 ;
+ fGains = 0 ;
+ fCalibrVersion= "v1" ;
+ fCalibrRun = 1 ;
+ fPatterns = 0 ;
+
TString clusterizerName( GetName()) ;
if (clusterizerName.IsNull() )
clusterizerName = "Default" ;
clusterizerName.Append(Version()) ;
SetName(clusterizerName) ;
fRecPointsInRun = 0 ;
- fCalibrationDB = 0 ;
}
Int_t index ;
//Evaluate position, dispersion and other RecPoint properties...
- for(index = 0; index < emcRecPoints->GetEntries(); index++)
- ((AliPHOSEmcRecPoint *)emcRecPoints->At(index))->EvalAll(fW0,digits) ;
-
+ for(index = 0; index < emcRecPoints->GetEntries(); index++){
+ AliPHOSEmcRecPoint * emcrp = static_cast<AliPHOSEmcRecPoint *>(emcRecPoints->At(index)) ;
+ emcrp->Purify(fPurifyThreshold) ;
+ if(emcrp->GetMultiplicity())
+ emcrp->EvalAll(fW0,digits) ;
+ else
+ emcRecPoints->Remove(emcrp) ;
+ }
+
+ emcRecPoints->Compress() ;
emcRecPoints->Sort() ;
for(index = 0; index < emcRecPoints->GetEntries(); index++)
((AliPHOSEmcRecPoint *)emcRecPoints->At(index))->SetIndexInList(index) ;
TClonesArray * digits = gime->Digits() ;
if ( !digits ) {
- Error("MakeClusters", "Digits with name %s not found !") ;
- abort() ;
+ Fatal("MakeClusters", "Digits with name %s not found !", BranchName().Data() ) ;
}
TClonesArray * digitsC = (TClonesArray*)digits->Clone() ;
TArrayI clusterdigitslist(1500) ;
Int_t index ;
+ // cout << "digit " << digit << " " << Calibrate(digit->GetAmp(),digit->GetId()) << endl ;
if (( IsInEmc (digit) && Calibrate(digit->GetAmp(),digit->GetId()) > fEmcClusteringThreshold ) ||
( IsInCpv (digit) && Calibrate(digit->GetAmp(),digit->GetId()) > fCpvClusteringThreshold ) ) {
if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
- emcRecPoints->Remove(emcRecPoint);
- emcRecPoints->Compress() ;
- index-- ;
- fNumberOfEmcClusters -- ;
- numberofNotUnfolded-- ;
+ if(emcRecPoint->GetNExMax()!=-1){ //Do not remove clusters which failed to unfold
+ emcRecPoints->Remove(emcRecPoint);
+ emcRecPoints->Compress() ;
+ index-- ;
+ fNumberOfEmcClusters -- ;
+ numberofNotUnfolded-- ;
+ }
}
+ else
+ emcRecPoint->SetNExMax(1) ; //Only one local maximum
delete[] maxAt ;
delete[] maxAtEnergy ;
Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
if( !rv ) {
- // Fit failed, return and remove cluster
+ // Fit failed, return
+ iniEmc->SetNExMax(-1) ;
delete[] fitparameters ;
return ;
}
(*emcRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
emcRP = (AliPHOSEmcRecPoint *) emcRecPoints->At(fNumberOfEmcClusters);
fNumberOfEmcClusters++ ;
+ emcRP->SetNExMax((Int_t)nPar/3) ;
}
else{//create new entries in fCpvRecPoints
if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
efit += x[iParam] * ShowerShape(distance) ;
iParam++ ;
}
- Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
+ Double_t sum=0 ;
+ if(emcEnergies[iDigit] > 0)
+ sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ;
+ // Here we assume, that sigma = sqrt(E)
iParam = 0 ;
while(iParam < nPar ){
Double_t xpar = x[iParam] ;
efit += epar * ShowerShape(distance) ;
}
- fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
+ if(emcEnergies[iDigit] > 0) //for the case if rounding error
+ fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/
+ emcEnergies[iDigit] ;
// Here we assume, that sigma = sqrt(E)
}
TObjArray * emcRecPoints = AliPHOSGetter::GetInstance()->EmcRecPoints(BranchName()) ;
TObjArray * cpvRecPoints = AliPHOSGetter::GetInstance()->CpvRecPoints(BranchName()) ;
- TString message ;
- message = "\nevent " ;
- message += gAlice->GetEvNumber() ;
- message += "\n Found " ;
- message += emcRecPoints->GetEntriesFast() ;
- message += "EMC RecPoints and " ;
- message += cpvRecPoints->GetEntriesFast() ;
- message += " CPV RecPoints \n" ;
-
- fRecPointsInRun += emcRecPoints->GetEntriesFast() ;
- fRecPointsInRun += cpvRecPoints->GetEntriesFast() ;
+ if(strstr(option,"deb")){
+ printf("event %d, found %d EmcRecPoints and %d CPV RecPoints \n",
+ gAlice->GetEvNumber(),
+ emcRecPoints->GetEntriesFast(),
+ cpvRecPoints->GetEntriesFast()) ;
+ }
+
if(strstr(option,"all")) {
- message += "EMC clusters " ;
- message += " Index Ene(MeV) Multi Module X Y Z Lambda 1 Lambda 2 # of prim Primaries list " ;
+ printf("\nEMC clusters\n" );
+ printf("Index Ene(MeV) Multi Module X Y Z Lambda 1 Lambda 2 # of prim Primaries list\n" );
Int_t index ;
for (index = 0 ; index < emcRecPoints->GetEntries() ; index++) {
AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )emcRecPoints->At(index) ;
rp->GetElipsAxis(lambda);
Int_t * primaries;
Int_t nprimaries;
- message += " " ;
primaries = rp->GetPrimaries(nprimaries);
- message += "\n" ;
- message += rp->GetIndexInList() ;
- message += " " ;
- message += rp->GetEnergy();
- message += " " ;
- message += rp->GetMultiplicity() ;
- message += " " ;
- message += rp->GetPHOSMod() ;
- message += " " ;
- message += locpos.X() ;
- message += " " ;
- message += locpos.Y() ;
- message += " " ;
- message += locpos.Z() ;
- message += " " ;
- message += lambda[0] ;
- message += " " ;
- message += lambda[1] ;
- message += " " ;
- message += nprimaries ;
- message += " : " ;
+ printf("\n%6d %8.2f %3d %2d %4.1f %4.1f %4.1f %4f %4f %2d : ",
+ rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetPHOSMod(),
+ locpos.X(), locpos.Y(), locpos.Z(), lambda[0], lambda[1], nprimaries) ;
for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
- message += primaries[iprimary] ;
- message += " " ;
+ printf("%d ", primaries[iprimary] ) ;
}
- }
-
+ }
//Now plot CPV recPoints
- message += "\nCPV clusters Index Multi Module X Y Z # of prim Primaries list \n" ;
+ printf("\n CPV RecPoints \n") ;
+ printf("Index Ene(MeV) Module X Y Z # of prim Primaries list\n") ;
for (index = 0 ; index < cpvRecPoints->GetEntries() ; index++) {
- AliPHOSRecPoint * rp = (AliPHOSRecPoint * )cpvRecPoints->At(index) ;
-
- TVector3 locpos;
- rp->GetLocalPosition(locpos);
-
- Int_t * primaries;
- Int_t nprimaries ;
- message += "\n" ;
- message += rp->GetIndexInList() ;
- message += " " ;
- message += rp->GetPHOSMod();
- message += " " ;
- message += locpos.X() ;
- message += " " ;
- message += locpos.Y() ;
- message += " " ;
- message += locpos.Z() ;
- message += " " ;
- message += nprimaries ;
- message += " : " ;
-
- primaries = rp->GetPrimaries(nprimaries);
- for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
- message += primaries[iprimary] ;
- message += " " ;
+ AliPHOSRecPoint * rp = (AliPHOSRecPoint * )cpvRecPoints->At(index) ;
+
+ TVector3 locpos;
+ rp->GetLocalPosition(locpos);
+
+ Int_t * primaries;
+ Int_t nprimaries ;
+ primaries = rp->GetPrimaries(nprimaries);
+ printf("\n%6d %8.2f %2d %4.1f %4.1f %4.1f %2d : ",
+ rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
+ locpos.X(), locpos.Y(), locpos.Z(), nprimaries) ;
+ primaries = rp->GetPrimaries(nprimaries);
+ for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
+ printf("%d ", primaries[iprimary] ) ;
+ }
}
- }
}
- Info("Print", message.Data() ) ;
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff