// This TTask is normally called from Reconstructioner, but can as well be used in
// standalone mode.
// Use Case:
-// root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1("galice.root")
+// root [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1("galice.root", "recpointsname", "digitsname")
// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
-// //reads gAlice from header file "..."
+// // reads gAlice from header file "galice.root", uses digits stored in the branch names "digitsname" (default = "Default")
+// // and saves recpoints in branch named "recpointsname" (default = "digitsname")
// root [1] cl->ExecuteTask()
// //finds RecPoints in all events stored in galice.root
// root [2] cl->SetDigitsBranch("digits2")
// --- Standard library ---
-#include <iostream.h>
-#include <iomanip.h>
-
// --- AliRoot header files ---
-
#include "AliPHOSClusterizerv1.h"
#include "AliPHOSCpvRecPoint.h"
#include "AliPHOSDigit.h"
#include "AliPHOSDigitizer.h"
#include "AliPHOSEmcRecPoint.h"
#include "AliPHOS.h"
-#include "AliPHOSPpsdRecPoint.h"
#include "AliPHOSGetter.h"
#include "AliRun.h"
+#include "AliPHOSCalibrManager.h"
+#include "AliPHOSCalibrationData.h"
ClassImp(AliPHOSClusterizerv1)
{
// default ctor (to be used mainly by Streamer)
- fNumberOfCpvClusters = 0 ;
- fNumberOfEmcClusters = 0 ;
-
- fCpvClusteringThreshold = 0.0;
- fEmcClusteringThreshold = 0.2;
- fPpsdClusteringThreshold = 0.0000002 ;
-
- fEmcLocMaxCut = 0.03 ;
- fCpvLocMaxCut = 0.03 ;
-
- fW0 = 4.5 ;
- fW0CPV = 4.0 ;
-
- fHeaderFileName = "" ;
- fDigitsBranchTitle = "" ;
-
+ InitParameters() ;
+ fDefaultInit = kTRUE ;
}
//____________________________________________________________________________
-AliPHOSClusterizerv1::AliPHOSClusterizerv1(const char* headerFile,const char* name)
-:AliPHOSClusterizer(headerFile, name)
+AliPHOSClusterizerv1::AliPHOSClusterizerv1(const char* headerFile,const char* name, const Bool_t toSplit)
+:AliPHOSClusterizer(headerFile, name, toSplit)
{
// ctor with the indication of the file where header Tree and digits Tree are stored
+ InitParameters() ;
+ Init() ;
+ fDefaultInit = kFALSE ;
- fNumberOfCpvClusters = 0 ;
- fNumberOfEmcClusters = 0 ;
-
- fCpvClusteringThreshold = 0.0;
- fEmcClusteringThreshold = 0.2;
- fPpsdClusteringThreshold = 0.0000002 ;
-
- fEmcLocMaxCut = 0.03 ;
- fCpvLocMaxCut = 0.03 ;
-
- fW0 = 4.5 ;
- fW0CPV = 4.0 ;
-
- fToUnfold = kTRUE ;
-
- fHeaderFileName = GetTitle() ;
- fDigitsBranchTitle = GetName() ;
+}
- TString tempo(GetName()) ;
- tempo.Append(Version()) ;
- SetName(tempo.Data()) ;
-
- Init() ;
+//____________________________________________________________________________
+ AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
+{
+ // dtor
+
+ delete fPedestals ;
+ delete fGains ;
+
+ fSplitFile = 0 ;
}
+//____________________________________________________________________________
+const TString AliPHOSClusterizerv1::BranchName() const
+{
+ // Returns branch name of the AliPHOSClusterizer
+
+ TString branchName(GetName() ) ;
+ branchName.Remove(branchName.Index(Version())-1) ;
+ return branchName ;
+}
+
+//____________________________________________________________________________
+Float_t AliPHOSClusterizerv1::Calibrate(Int_t amp, Int_t absId) const
+{
+ // Returns energy value in the cell absId with the digit amplitude amp
+
+ 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 ;
+ else //calibrate as CPV
+ return fADCpedestalCpv+ amp*fADCchanelCpv ;
+ }
+}
+
//____________________________________________________________________________
void AliPHOSClusterizerv1::Exec(Option_t * option)
{
// Steering method
-
if( strcmp(GetName(), "")== 0 )
Init() ;
if(strstr(option,"tim"))
- gBenchmark->Start("PHOSClusterizerv1");
+ gBenchmark->Start("PHOSClusterizer");
if(strstr(option,"print"))
Print("") ;
- //check, if the branch with name of this" already exits?
- TObjArray * lob = (TObjArray*)gAlice->TreeR()->GetListOfBranches() ;
- TIter next(lob) ;
- TBranch * branch = 0 ;
- Bool_t phosemcfound = kFALSE, phoscpvfound = kFALSE, clusterizerfound = kFALSE ;
-
- TString taskName(GetName()) ;
- taskName.ReplaceAll(Version(), "") ;
+ AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
+ if(gime->BranchExists("RecPoints"))
+ return ;
- while ( (branch = (TBranch*)next()) && (!phosemcfound || !phoscpvfound || !clusterizerfound) ) {
- if ( (strcmp(branch->GetName(), "PHOSEmcRP")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) )
- phosemcfound = kTRUE ;
-
- else if ( (strcmp(branch->GetName(), "PHOSCpvRP")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) )
- phoscpvfound = kTRUE ;
-
- else if ( (strcmp(branch->GetName(), "AliPHOSClusterizer")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) )
- clusterizerfound = kTRUE ;
+ //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) ;
}
- if ( phoscpvfound || phosemcfound || clusterizerfound ) {
- cerr << "WARNING: AliPHOSClusterizer::Exec -> Emc(Cpv)RecPoints and/or Clusterizer branch with name "
- << taskName.Data() << " already exits" << endl ;
- return ;
- }
-
- Int_t nevents = (Int_t) gAlice->TreeE()->GetEntries() ;
+ 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() ;
- if(!ReadDigits(ievent)) //reads digits for event fEvent
- return;
-
+ fNumberOfEmcClusters = fNumberOfCpvClusters = 0 ;
+
MakeClusters() ;
- if(fToUnfold)
+ if(fToUnfold)
MakeUnfolding() ;
WriteRecPoints(ievent) ;
- if(strstr(option,"deb"))
+ if(strstr(option,"deb"))
PrintRecPoints(option) ;
- }
+
+ //increment the total number of digits per run
+ fRecPointsInRun += gime->EmcRecPoints()->GetEntriesFast() ;
+ fRecPointsInRun += gime->CpvRecPoints()->GetEntriesFast() ;
+ }
if(strstr(option,"tim")){
gBenchmark->Stop("PHOSClusterizer");
- cout << "AliPHOSClusterizer:" << endl ;
- cout << " took " << gBenchmark->GetCpuTime("PHOSClusterizer") << " seconds for Clusterizing "
- << gBenchmark->GetCpuTime("PHOSClusterizer")/nevents << " seconds per event " << endl ;
- cout << endl ;
+ Info("Exec", " took %f seconds for Clusterizing %f seconds per event \n",
+ gBenchmark->GetCpuTime("PHOSClusterizer"),
+ gBenchmark->GetCpuTime("PHOSClusterizer")/nevents ) ;
}
-
}
//____________________________________________________________________________
-Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, int * maxAt, Float_t * maxAtEnergy,
+Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** maxAt, Float_t * maxAtEnergy,
Int_t nPar, Float_t * fitparameters) const
{
// Calls TMinuit to fit the energy distribution of a cluster with several maxima
const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
for(iDigit = 0; iDigit < nDigits; iDigit++){
- digit = (AliPHOSDigit *) maxAt[iDigit];
+ digit = maxAt[iDigit];
Int_t relid[4] ;
Float_t x = 0.;
gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ;
index++ ;
if(ierflg != 0){
- cout << "PHOS Unfolding> Unable to set initial value for fit procedure : x = " << x << endl ;
+ Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : x = %f\n", x ) ;
return kFALSE;
}
gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ;
index++ ;
if(ierflg != 0){
- cout << "PHOS Unfolding> Unable to set initial value for fit procedure : z = " << z << endl ;
+ Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : z =%f\n", z ) ;
return kFALSE;
}
gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ;
index++ ;
if(ierflg != 0){
- cout << "PHOS Unfolding> Unable to set initial value for fit procedure : energy = " << energy << endl ;
+ Warning("FindFit", "PHOS Unfolding unable to set initial value for fit procedure : energy = %f\n", energy ) ;
return kFALSE;
}
}
gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize
if(ierflg == 4){ // Minimum not found
- cout << "PHOS Unfolding> Fit not converged, cluster abandoned "<< endl ;
+ Warning("FindFit", "PHOS Unfolding fit not converged, cluster abandoned\n" );
return kFALSE ;
}
for(index = 0; index < nPar; index++){
}
+//____________________________________________________________________________
+void AliPHOSClusterizerv1::GetCalibrationParameters()
+{
+ // Get calibration parameters from AliPHOSDigitizer
+
+ AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
+
+ const TTask * task = gime->Digitizer(BranchName()) ;
+ if(strcmp(task->IsA()->GetName(),"AliPHOSDigitizer")==0){
+ const AliPHOSDigitizer * dig = static_cast<const AliPHOSDigitizer *>(task) ;
+
+ fADCchanelEmc = dig->GetEMCchannel() ;
+ fADCpedestalEmc = dig->GetEMCpedestal();
+
+ fADCchanelCpv = dig->GetCPVchannel() ;
+ fADCpedestalCpv = dig->GetCPVpedestal() ;
+ }
+ else{
+ ;
+// fCalibrationDB = gime->CalibrationDB();
+ }
+}
+
//____________________________________________________________________________
void AliPHOSClusterizerv1::Init()
{
if ( strcmp(GetTitle(), "") == 0 )
SetTitle("galice.root") ;
- TString taskName(GetName()) ;
- taskName.ReplaceAll(Version(), "") ;
+ TString branchname = GetName() ;
+ branchname.Remove(branchname.Index(Version())-1) ;
- AliPHOSGetter * gime = AliPHOSGetter::GetInstance(GetTitle(), taskName.Data()) ;
+ AliPHOSGetter * gime = AliPHOSGetter::GetInstance(GetTitle(),branchname.Data(), fToSplit ) ;
if ( gime == 0 ) {
- cerr << "ERROR: AliPHOSClusterizerv1::Init -> Could not obtain the Getter object !" << endl ;
+ Error("Init", "Could not obtain the Getter object !") ;
return ;
}
+
+ fSplitFile = 0 ;
+ if(fToSplit){
+ // construct the name of the file as /path/EMCAL.SDigits.root
+ //First - extract full path if necessary
+ TString fileName(GetTitle()) ;
+ Ssiz_t islash = fileName.Last('/') ;
+ if(islash<fileName.Length())
+ fileName.Remove(islash+1,fileName.Length()) ;
+ else
+ fileName="" ;
+ // Next - append the file name
+ fileName+="PHOS.RecData." ;
+ if((strcmp(branchname.Data(),"Default")!=0)&&(strcmp(branchname.Data(),"")!=0)){
+ fileName+=branchname ;
+ fileName+="." ;
+ }
+ fileName+="root" ;
+ // Finally - check if the file already opened or open the file
+ fSplitFile = static_cast<TFile*>(gROOT->GetFile(fileName.Data()));
+ if(!fSplitFile)
+ fSplitFile = TFile::Open(fileName.Data(),"update") ;
+ }
+
+
+ const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
+ fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ;
+
if(!gMinuit)
gMinuit = new TMinuit(100) ;
- //add Task to //YSAlice/tasks/Reconstructioner/PHOS
- TTask * aliceRe = (TTask*)gROOT->FindObjectAny("YSAlice/tasks/Reconstructioner") ;
- TTask * phosRe = (TTask*)aliceRe->GetListOfTasks()->FindObject("PHOS") ;
- phosRe->Add(this) ;
- // create a folder on the white board //YSAlice/WhiteBoard/RecPoints/PHOS/recpointsName
- gime->Post(GetTitle(), "R", taskName.Data() ) ;
+ gime->PostClusterizer(this) ;
+ gime->PostRecPoints(branchname) ;
}
+//____________________________________________________________________________
+void AliPHOSClusterizerv1::InitParameters()
+{
+ // Set initial parameters for clusterization
+
+ fNumberOfCpvClusters = 0 ;
+ fNumberOfEmcClusters = 0 ;
+
+ fCpvClusteringThreshold = 0.0;
+ fEmcClusteringThreshold = 0.2;
+
+ fEmcLocMaxCut = 0.03 ;
+ fCpvLocMaxCut = 0.03 ;
+
+ fW0 = 4.5 ;
+ fW0CPV = 4.0 ;
+
+ 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(":") ;
+ clusterizerName.Append(Version()) ;
+ SetName(clusterizerName) ;
+ fRecPointsInRun = 0 ;
+
+
+}
+
//____________________________________________________________________________
Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
{
Int_t relid2[4] ;
geom->AbsToRelNumbering(d2->GetId(), relid2) ;
- if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module and the same PPSD Module
+ if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
Int_t coldiff = TMath::Abs( relid1[3] - relid2[3] ) ;
if (( coldiff <= 1 ) && ( rowdiff <= 1 )){
+ if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate))
rv = 1 ;
}
else {
}
else {
- if( (relid1[0] < relid2[0]) || (relid1[1] < relid2[1]) )
+ if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) )
rv=2 ;
}
- //Do NOT clusterize upper PPSD
- if( IsInPpsd(d1) && IsInPpsd(d2) &&
- relid1[1] > 0 &&
- relid1[1] < geom->GetNumberOfPadsPhi()*geom->GetNumberOfPadsPhi() ) rv = 2 ;
-
return rv ;
}
Bool_t rv = kFALSE ;
const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ;
- Int_t relid[4] ;
- geom->AbsToRelNumbering(digit->GetId(), relid) ;
+ Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
- if ( relid[1] == 0 ) rv = kTRUE;
-
- return rv ;
-}
-
-//____________________________________________________________________________
-Bool_t AliPHOSClusterizerv1::IsInPpsd(AliPHOSDigit * digit) const
-{
- // Tells if (true) or not (false) the digit is in a PHOS-PPSD module
-
- Bool_t rv = kFALSE ;
- const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ;
-
- Int_t relid[4] ;
- geom->AbsToRelNumbering(digit->GetId(), relid) ;
-
- if ( relid[1] > 0 && relid[0] > geom->GetNCPVModules() ) rv = kTRUE;
+ if(digit->GetId() <= nEMC ) rv = kTRUE;
return rv ;
}
Bool_t rv = kFALSE ;
const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ;
- Int_t relid[4] ;
- geom->AbsToRelNumbering(digit->GetId(), relid) ;
+ Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
- if ( relid[1] > 0 && relid[0] <= geom->GetNCPVModules() ) rv = kTRUE;
+ if(digit->GetId() > nEMC ) rv = kTRUE;
return rv ;
}
-//____________________________________________________________________________
-Bool_t AliPHOSClusterizerv1::ReadDigits(Int_t event)
- {
- // reads digitis with specified title from TreeD
-
- fNumberOfEmcClusters = 0 ;
- fNumberOfCpvClusters = 0 ;
-
- // Get Digits Tree header from file
- gAlice->GetEvent(event) ;
- gAlice->SetEvent(event) ;
-
- if ( gAlice->TreeD()==0) {
- cerr << "ERROR: AliPHOSClusterizerv1::ReadDigits There is no Digit Tree" << endl;
- return kFALSE;
- }
-
- //set address of the Digits and Digitizer
- TBranch * digitsBranch = 0;
- TBranch * digitizerBranch = 0;
- TObjArray * lob = (TObjArray*)gAlice->TreeD()->GetListOfBranches() ;
- TIter next(lob) ;
- TBranch * branch = 0 ;
- Bool_t phosfound = kFALSE, digitizerfound = kFALSE ;
-
- TString taskName(GetName()) ;
- taskName.ReplaceAll(Version(), "") ;
-
- while ( (branch = (TBranch*)next()) && (!phosfound || !digitizerfound) ) {
- if ( (strcmp(branch->GetName(), "PHOS")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) ) {
- phosfound = kTRUE ;
- digitsBranch = branch ;
- }
-
- else if ( (strcmp(branch->GetName(), "AliPHOSDigitizer")==0) && (strcmp(branch->GetTitle(), taskName.Data())==0) ) {
- digitizerfound = kTRUE ;
- digitizerBranch = branch ;
- }
- }
- if ( !phosfound || !digitizerfound ) {
- cerr << "WARNING: AliPHOSClusterizerv1::ReadDigits -> Digits and/or Digitizer branch with name " << taskName.Data()
- << " not found" << endl ;
- return kFALSE ;
- }
-
- AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
-
- TClonesArray * digits = gime->Digits() ;
- digits->Clear() ;
- digitsBranch->SetAddress(&digits) ;
-
- AliPHOSDigitizer * digitizer = gime->Digitizer() ;
- digitizerBranch->SetAddress(&digitizer) ;
-
- digitsBranch->GetEntry(0) ;
- digitizerBranch->GetEntry(0) ;
-
- fPedestal = digitizer->GetPedestal() ;
- fSlope = digitizer->GetSlope() ;
-
- return kTRUE ;
-}
//____________________________________________________________________________
void AliPHOSClusterizerv1::WriteRecPoints(Int_t event)
// Creates new branches with given title
// fills and writes into TreeR.
+
AliPHOSGetter *gime = AliPHOSGetter::GetInstance() ;
TObjArray * emcRecPoints = gime->EmcRecPoints() ;
TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
TClonesArray * digits = gime->Digits() ;
+ TTree * treeR ;
+
+ if(fToSplit){
+ if(!fSplitFile)
+ return ;
+ fSplitFile->cd() ;
+ TString name("TreeR") ;
+ name += event ;
+ treeR = dynamic_cast<TTree*>(fSplitFile->Get(name));
+ }
+ else{
+ treeR = gAlice->TreeR();
+ }
+
+ if(!treeR){
+ if(fSplitFile)
+ gAlice->MakeTree("R", fSplitFile);
+ else
+ gAlice->MakeTree("R", gROOT->GetFile(GetTitle()));
+ treeR = gAlice->TreeR() ;
+ }
Int_t index ;
- //Evaluate poisition, dispersion and other RecPoint properties...
- for(index = 0; index < emcRecPoints->GetEntries(); index++)
- ((AliPHOSEmcRecPoint *)emcRecPoints->At(index))->EvalAll(fW0,digits) ;
+ //Evaluate position, dispersion and other RecPoint properties...
+ 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) ;
-
+
emcRecPoints->Expand(emcRecPoints->GetEntriesFast()) ;
-
+
//Now the same for CPV
for(index = 0; index < cpvRecPoints->GetEntries(); index++)
((AliPHOSRecPoint *)cpvRecPoints->At(index))->EvalAll(fW0CPV,digits) ;
-
+
cpvRecPoints->Sort() ;
-
+
for(index = 0; index < cpvRecPoints->GetEntries(); index++)
((AliPHOSRecPoint *)cpvRecPoints->At(index))->SetIndexInList(index) ;
-
- cpvRecPoints->Expand(cpvRecPoints->GetEntriesFast()) ;
-
- gAlice->GetEvent(event) ;
- if(gAlice->TreeR()==0)
- gAlice->MakeTree("R") ;
- //Make branches in TreeR for RecPoints and Clusterizer
- char * filename = 0;
- if(gSystem->Getenv("CONFIG_SPLIT_FILE")!=0){ //generating file name
- filename = new char[strlen(gAlice->GetBaseFile())+20] ;
- sprintf(filename,"%s/PHOS.Reco.root",gAlice->GetBaseFile()) ;
- }
-
- //Make new branches
- TDirectory *cwd = gDirectory;
+ cpvRecPoints->Expand(cpvRecPoints->GetEntriesFast()) ;
- TString branchName(GetName()) ;
- branchName.ReplaceAll(Version(), "") ;
-
- //First EMC
Int_t bufferSize = 32000 ;
Int_t splitlevel = 0 ;
- TBranch * emcBranch = gAlice->TreeR()->Branch("PHOSEmcRP","TObjArray",&emcRecPoints,bufferSize,splitlevel);
- emcBranch->SetTitle(branchName);
- if (filename) {
- emcBranch->SetFile(filename);
- TIter next( emcBranch->GetListOfBranches());
- TBranch * sb ;
- while ((sb=(TBranch*)next())) {
- sb->SetFile(filename);
- }
-
- cwd->cd();
- }
+
+ //First EMC
+ TBranch * emcBranch = treeR->Branch("PHOSEmcRP","TObjArray",&emcRecPoints,bufferSize,splitlevel);
+ emcBranch->SetTitle(BranchName());
//Now CPV branch
- TBranch * cpvBranch = gAlice->TreeR()->Branch("PHOSCpvRP","TObjArray",&cpvRecPoints,bufferSize,splitlevel);
- cpvBranch->SetTitle(branchName);
- if (filename) {
- cpvBranch->SetFile(filename);
- TIter next( cpvBranch->GetListOfBranches());
- TBranch * sb;
- while ((sb=(TBranch*)next())) {
- sb->SetFile(filename);
- }
- cwd->cd();
- }
+ TBranch * cpvBranch = treeR->Branch("PHOSCpvRP","TObjArray",&cpvRecPoints,bufferSize,splitlevel);
+ cpvBranch->SetTitle(BranchName());
//And Finally clusterizer branch
- AliPHOSClusterizerv1 * cl = (AliPHOSClusterizerv1*)gime->Clusterizer(GetName()) ;
- TBranch * clusterizerBranch = gAlice->TreeR()->Branch("AliPHOSClusterizer","AliPHOSClusterizerv1",
+ AliPHOSClusterizerv1 * cl = this ;
+ TBranch * clusterizerBranch = treeR->Branch("AliPHOSClusterizer","AliPHOSClusterizerv1",
&cl,bufferSize,splitlevel);
- clusterizerBranch->SetTitle(branchName);
- if (filename) {
- clusterizerBranch->SetFile(filename);
- TIter next( clusterizerBranch->GetListOfBranches());
- TBranch * sb ;
- while ((sb=(TBranch*)next())) {
- sb->SetFile(filename);
- }
- cwd->cd();
- }
-
- emcBranch->Fill() ;
- cpvBranch->Fill() ;
+ clusterizerBranch->SetTitle(BranchName());
+
+ emcBranch ->Fill() ;
+ cpvBranch ->Fill() ;
clusterizerBranch->Fill() ;
- gAlice->TreeR()->Write(0,kOverwrite) ;
-
+ treeR->AutoSave() ;
+ if(gAlice->TreeR()!=treeR)
+ treeR->Delete();
}
//____________________________________________________________________________
AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
- TObjArray * emcRecPoints = gime->EmcRecPoints() ;
- TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
- emcRecPoints->Clear() ;
- cpvRecPoints->Clear() ;
+ TObjArray * emcRecPoints = gime->EmcRecPoints(BranchName()) ;
+ TObjArray * cpvRecPoints = gime->CpvRecPoints(BranchName()) ;
+ emcRecPoints->Delete() ;
+ cpvRecPoints->Delete() ;
TClonesArray * digits = gime->Digits() ;
+ if ( !digits ) {
+ 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()) > fEmcClusteringThreshold ) ||
- ( IsInPpsd(digit) && Calibrate(digit->GetAmp()) > fPpsdClusteringThreshold ) ||
- ( IsInCpv (digit) && Calibrate(digit->GetAmp()) > fCpvClusteringThreshold ) ) {
-
+ if (( IsInEmc (digit) && Calibrate(digit->GetAmp(),digit->GetId()) > fEmcClusteringThreshold ) ||
+ ( IsInCpv (digit) && Calibrate(digit->GetAmp(),digit->GetId()) > fCpvClusteringThreshold ) ) {
Int_t iDigitInCluster = 0 ;
if ( IsInEmc(digit) ) {
if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
emcRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
- emcRecPoints->AddAt(new AliPHOSEmcRecPoint(), fNumberOfEmcClusters) ;
+ emcRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
clu = (AliPHOSEmcRecPoint *) emcRecPoints->At(fNumberOfEmcClusters) ;
fNumberOfEmcClusters++ ;
- clu->AddDigit(*digit, Calibrate(digit->GetAmp())) ;
+ clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->GetId())) ;
clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
iDigitInCluster++ ;
digitsC->Remove(digit) ;
} else {
- // start a new PPSD/CPV cluster
+ // start a new CPV cluster
if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
cpvRecPoints->Expand(2*fNumberOfCpvClusters+1);
- if(IsInPpsd(digit))
- cpvRecPoints->AddAt(new AliPHOSPpsdRecPoint(),fNumberOfCpvClusters) ;
- else
- cpvRecPoints->AddAt(new AliPHOSCpvRecPoint(), fNumberOfCpvClusters) ;
+ cpvRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
- clu = (AliPHOSPpsdRecPoint *) cpvRecPoints->At(fNumberOfCpvClusters) ;
+ clu = (AliPHOSCpvRecPoint *) cpvRecPoints->At(fNumberOfCpvClusters) ;
fNumberOfCpvClusters++ ;
- clu->AddDigit(*digit, Calibrate(digit->GetAmp()) ) ;
+ clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->GetId()) ) ;
clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
iDigitInCluster++ ;
digitsC->Remove(digit) ;
case 0 : // not a neighbour
break ;
case 1 : // are neighbours
- clu->AddDigit(*digitN, Calibrate( digitN->GetAmp() ) ) ;
+ clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), digitN->GetId() ) ) ;
clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ;
iDigitInCluster++ ;
digitsC->Remove(digitN) ;
void AliPHOSClusterizerv1::MakeUnfolding()
{
// Unfolds clusters using the shape of an ElectroMagnetic shower
- // Performs unfolding of all EMC/CPV but NOT ppsd clusters
+ // Performs unfolding of all EMC/CPV clusters
AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
+
const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
- TObjArray * emcRecPoints = gime->EmcRecPoints() ;
- TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
+ TObjArray * emcRecPoints = gime->EmcRecPoints(BranchName()) ;
+ TObjArray * cpvRecPoints = gime->CpvRecPoints(BranchName()) ;
TClonesArray * digits = gime->Digits() ;
// Unfold first EMC clusters
break ;
Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
- Int_t * maxAt = new Int_t[nMultipl] ;
+ AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
Float_t * maxAtEnergy = new Float_t[nMultipl] ;
Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,digits) ;
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 ;
// Unfold now CPV clusters
if(fNumberOfCpvClusters > 0){
- Int_t nModulesToUnfold = geom->GetNCPVModules() ;
+ Int_t nModulesToUnfold = geom->GetNModules() ;
Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
Int_t index ;
AliPHOSEmcRecPoint * emcRecPoint = (AliPHOSEmcRecPoint*) recPoint ;
Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
- Int_t * maxAt = new Int_t[nMultipl] ;
+ AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
Float_t * maxAtEnergy = new Float_t[nMultipl] ;
Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,digits) ;
//____________________________________________________________________________
void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
Int_t nMax,
- int * maxAt,
+ AliPHOSDigit ** maxAt,
Float_t * maxAtEnergy)
{
// Performs the unfolding of a cluster with nMax overlapping showers
AliPHOSGetter * gime = AliPHOSGetter::GetInstance() ;
+
const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
- TClonesArray * digits = gime->Digits() ;
+ const TClonesArray * digits = gime->Digits() ;
TObjArray * emcRecPoints = gime->EmcRecPoints() ;
TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
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 ;
}
if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
emcRecPoints->Expand(2*fNumberOfEmcClusters) ;
- (*emcRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint() ;
+ (*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())
cpvRecPoints->Expand(2*fNumberOfCpvClusters) ;
- (*cpvRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint() ;
+ (*cpvRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
emcRP = (AliPHOSEmcRecPoint *) cpvRecPoints->At(fNumberOfCpvClusters);
fNumberOfCpvClusters++ ;
}
// Calculates the Chi square for the cluster unfolding minimization
// Number of parameters, Gradient, Chi squared, parameters, what to do
-
TList * toMinuit = (TList*) gMinuit->GetObjectFit() ;
AliPHOSEmcRecPoint * emcRP = (AliPHOSEmcRecPoint*) toMinuit->At(0) ;
Float_t * emcEnergies = emcRP->GetEnergiesList() ;
- AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
-
+ const AliPHOSGeometry * geom = AliPHOSGetter::GetInstance()->PHOSGeometry() ;
fret = 0. ;
Int_t iparam ;
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)
}
{
// Print clusterizer parameters
- if( strcmp(GetName(), "") !=0 ){
-
+ TString message ;
+ TString taskName(GetName()) ;
+ taskName.ReplaceAll(Version(), "") ;
+
+ if( strcmp(GetName(), "") !=0 ) {
// Print parameters
-
- TString taskName(GetName()) ;
- taskName.ReplaceAll(Version(), "") ;
-
- cout << "---------------"<< taskName.Data() << " " << GetTitle()<< "-----------" << endl
- << "Clusterizing digits from the file: " << fHeaderFileName.Data() << endl
- << " Branch: " << fDigitsBranchTitle.Data() << endl
- << endl
- << " EMC Clustering threshold = " << fEmcClusteringThreshold << endl
- << " EMC Local Maximum cut = " << fEmcLocMaxCut << endl
- << " EMC Logarothmic weight = " << fW0 << endl
- << endl
- << " CPV Clustering threshold = " << fCpvClusteringThreshold << endl
- << " CPV Local Maximum cut = " << fCpvLocMaxCut << endl
- << " CPV Logarothmic weight = " << fW0CPV << endl
- << endl
- << " PPSD Clustering threshold = " << fPpsdClusteringThreshold << endl;
+ message = "\n--------------- %s %s -----------\n" ;
+ message += "Clusterizing digits from the file: %s\n" ;
+ message += " Branch: %s\n" ;
+ message += " EMC Clustering threshold = %f\n" ;
+ message += " EMC Local Maximum cut = %f\n" ;
+ message += " EMC Logarothmic weight = %f\n" ;
+ message += " CPV Clustering threshold = %f\n" ;
+ message += " CPV Local Maximum cut = %f\n" ;
+ message += " CPV Logarothmic weight = %f\n" ;
if(fToUnfold)
- cout << " Unfolding on " << endl ;
+ message += " Unfolding on\n" ;
else
- cout << " Unfolding off " << endl ;
+ message += " Unfolding off\n" ;
- cout << "------------------------------------------------------------------" <<endl ;
+ message += "------------------------------------------------------------------" ;
}
- else
- cout << " AliPHOSClusterizerv1 not initialized " << endl ;
+ else
+ message = " AliPHOSClusterizerv1 not initialized " ;
+
+ Info("Print", message.Data(),
+ taskName.Data(),
+ GetTitle(),
+ taskName.Data(),
+ GetName(),
+ fEmcClusteringThreshold,
+ fEmcLocMaxCut,
+ fW0,
+ fCpvClusteringThreshold,
+ fCpvLocMaxCut,
+ fW0CPV ) ;
}
+
+
//____________________________________________________________________________
void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
{
// Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
- TObjArray * emcRecPoints = AliPHOSGetter::GetInstance()->EmcRecPoints() ;
- TObjArray * cpvRecPoints = AliPHOSGetter::GetInstance()->CpvRecPoints() ;
-
- cout << "AliPHOSClusterizerv1: " << endl ;
- cout << " Found "<< emcRecPoints->GetEntriesFast() << " EMC Rec Points and "
- << cpvRecPoints->GetEntriesFast() << " CPV RecPoints" << endl ;
+ TObjArray * emcRecPoints = AliPHOSGetter::GetInstance()->EmcRecPoints(BranchName()) ;
+ TObjArray * cpvRecPoints = AliPHOSGetter::GetInstance()->CpvRecPoints(BranchName()) ;
- if(strstr(option,"all")) {
- cout << "EMC clusters " << endl ;
- cout << " Index "
- << " Ene(MeV) "
- << " Multi "
- << " Module "
- << " X "
- << " Y "
- << " Z "
- << " Lambda 1 "
- << " Lambda 2 "
- << " MaxEnergy "
- << " # of prim "
- << " Primaries list " << endl;
+ 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")) {
+ 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) ;
-
- cout << setw(6) << rp->GetIndexInList() << " ";
- cout << setw(6) << rp->GetEnergy() << " ";
- cout << setw(6) << rp->GetMultiplicity()<< " ";
- cout << setw(6) << rp->GetPHOSMod() << " ";
-
TVector3 locpos;
rp->GetLocalPosition(locpos);
- cout << setw(8) << locpos.X() << " ";
- cout << setw(8) << locpos.Y() << " ";
- cout << setw(8) << locpos.Z() << " ";
-
Float_t lambda[2];
rp->GetElipsAxis(lambda);
- cout << setw(10)<< lambda[0] << " ";
- cout << setw(10)<< lambda[1] << " ";
-
-
Int_t * primaries;
Int_t nprimaries;
primaries = rp->GetPrimaries(nprimaries);
- cout << setw(8) << primaries << " ";
-
- for (Int_t iprimary=0; iprimary<nprimaries; iprimary++)
- cout << setw(4) << primaries[iprimary] << " ";
- cout << endl;
- }
- cout << endl ;
-
- //Now plot CPV/PPSD recPoints
- cout << "EMC clusters " << endl ;
- cout << " Index "
- << " Multi "
- << " Module "
- << " Layer "
- << " X "
- << " Y "
- << " Z "
- << " # of prim "
- << " Primaries list " << endl;
-
+ 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++) {
+ printf("%d ", primaries[iprimary] ) ;
+ }
+ }
+ //Now plot CPV recPoints
+ 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) ;
- cout << setw(6) << rp->GetIndexInList() << " ";
- cout << setw(6) << rp->GetPHOSMod() << " ";
-
- if( (strcmp(rp->ClassName() , "AliPHOSPpsdRecPoint" )) == 0){
- AliPHOSPpsdRecPoint * ppsd = (AliPHOSPpsdRecPoint*) rp ;
- if(ppsd->GetUp())
- cout <<" CPV ";
- else
- cout <<" PPSD ";
+ 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] ) ;
+ }
}
- else
- cout <<" CPV ";
-
- TVector3 locpos;
- rp->GetLocalPosition(locpos);
- cout << setw(8) << locpos.X() << " ";
- cout << setw(8) << locpos.Y() << " ";
- cout << setw(8) << locpos.Z() << " ";
-
- Int_t * primaries;
- Int_t nprimaries;
- primaries = rp->GetPrimaries(nprimaries);
- cout << setw(8) << primaries << " ";
-
- for (Int_t iprimary=0; iprimary<nprimaries; iprimary++)
- cout << setw(4) << primaries[iprimary] << " ";
- cout << endl;
- }
-
-
- cout << "-------------------------------------------------"<<endl ;
}
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff