/* $Id$ */
-/* $Log:
- 1 October 2000. Yuri Kharlov:
- AreNeighbours()
- PPSD upper layer is considered if number of layers>1
-
- 18 October 2000. Yuri Kharlov:
- AliPHOSClusterizerv1()
- CPV clusterizing parameters added
-
- MakeClusters()
- After first PPSD digit remove EMC digits only once
-*/
+/* History of cvs commits:
+ *
+ * $Log$
+ * Revision 1.82 2005/09/02 15:43:13 kharlov
+ * Add comments in GetCalibrationParameters and Calibrate
+ *
+ * Revision 1.81 2005/09/02 14:32:07 kharlov
+ * Calibration of raw data
+ *
+ * Revision 1.80 2005/08/24 15:31:36 kharlov
+ * Setting raw digits flag
+ *
+ * Revision 1.79 2005/07/25 15:53:53 kharlov
+ * Read raw data
+ *
+ * Revision 1.78 2005/05/28 14:19:04 schutz
+ * Compilation warnings fixed by T.P.
+ *
+ */
+
//*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute)
//////////////////////////////////////////////////////////////////////////////
// Clusterization class. Performs clusterization (collects neighbouring active cells) and
// 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")
// --- ROOT system ---
-#include "TROOT.h"
-#include "TFile.h"
-#include "TFolder.h"
#include "TMath.h"
#include "TMinuit.h"
#include "TTree.h"
-#include "TSystem.h"
#include "TBenchmark.h"
// --- Standard library ---
-#include <iostream.h>
-#include <iomanip.h>
-
// --- AliRoot header files ---
-
+#include "AliLog.h"
+#include "AliPHOSGetter.h"
+#include "AliPHOSGeometry.h"
#include "AliPHOSClusterizerv1.h"
+#include "AliPHOSEmcRecPoint.h"
#include "AliPHOSCpvRecPoint.h"
#include "AliPHOSDigit.h"
#include "AliPHOSDigitizer.h"
-#include "AliPHOSEmcRecPoint.h"
-#include "AliPHOS.h"
-#include "AliPHOSPpsdRecPoint.h"
-#include "AliRun.h"
+#include "AliPHOSCalibrationDB.h"
ClassImp(AliPHOSClusterizerv1)
//____________________________________________________________________________
- AliPHOSClusterizerv1::AliPHOSClusterizerv1():AliPHOSClusterizer()
+ AliPHOSClusterizerv1::AliPHOSClusterizerv1() : AliPHOSClusterizer()
{
// default ctor (to be used mainly by Streamer)
- SetName("AliPHOSClusterizer");
- SetTitle("Version 1") ;
-
- 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 ;
-
- fGeom = 0 ;
-
- fDigits = 0 ;
- fDigitizer = 0 ;
- fEmcRecPoints = 0 ;
- fCpvRecPoints = 0 ;
-
- fIsInitialized = kFALSE ;
+ InitParameters() ;
+ fDefaultInit = kTRUE ;
}
+
//____________________________________________________________________________
-AliPHOSClusterizerv1::AliPHOSClusterizerv1(const char* headerFile,const char* digitsFile):AliPHOSClusterizer()
+AliPHOSClusterizerv1::AliPHOSClusterizerv1(const TString alirunFileName, const TString eventFolderName)
+:AliPHOSClusterizer(alirunFileName, eventFolderName)
{
// ctor with the indication of the file where header Tree and digits Tree are stored
- SetName("AliPHOSClusterizer");
- SetTitle("Version 1") ;
-
- 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 = headerFile ;
- fDigitsBranchTitle = digitsFile ;
-
- TFile * file = (TFile*) gROOT->GetFile(fHeaderFileName.Data() ) ;
-
- if(file == 0){
- if(fHeaderFileName.Contains("rfio")) // if we read file using HPSS
- file = TFile::Open(fHeaderFileName.Data(),"update") ;
- else
- file = new TFile(fHeaderFileName.Data(),"update") ;
- gAlice = (AliRun *) file->Get("gAlice") ;
- }
-
- AliPHOS * phos = (AliPHOS *) gAlice->GetDetector("PHOS") ;
- fGeom = AliPHOSGeometry::GetInstance(phos->GetGeometry()->GetName(),phos->GetGeometry()->GetTitle() );
-
- fDigits = new TClonesArray("AliPHOSDigit",10) ;
- fDigitizer = new AliPHOSDigitizer() ;
- fEmcRecPoints = new TObjArray(200) ;
- fCpvRecPoints = new TObjArray(200) ;
- fEmcRecPoints->SetOwner(); // This lets Clear() really detete rec.points in array
- fCpvRecPoints->SetOwner();
-
- if(!gMinuit) gMinuit = new TMinuit(100) ;
-
- // add Task to //root/Tasks folder
- TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ;
- TTask * phostasks = (TTask*)(roottasks->GetListOfTasks()->FindObject("PHOS"));
- phostasks->Add(this) ;
-
+ InitParameters() ;
+ Init() ;
+ fDefaultInit = kFALSE ;
+}
- fIsInitialized = kTRUE ;
+//____________________________________________________________________________
+ AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
+{
+ // dtor
}
//____________________________________________________________________________
-void AliPHOSClusterizerv1::Exec(Option_t * option)
-{
- // Steering function
+const TString AliPHOSClusterizerv1::BranchName() const
+{
+ return GetName();
+}
+
+//____________________________________________________________________________
+Float_t AliPHOSClusterizerv1::Calibrate(Int_t amp, Int_t absId)
+{
+ // Convert digitized amplitude into energy.
+ // Calibration parameters are taken from calibration data base for raw data,
+ // or from digitizer parameters for simulated data.
+
+ if(fCalibData){
+ Int_t relId[4];
+ AliPHOSGetter *gime = AliPHOSGetter::Instance();
+ gime->PHOSGeometry()->AbsToRelNumbering(absId,relId) ;
+ Int_t module = relId[0];
+ Int_t row = relId[2];
+ Int_t column = relId[3];
+ if(absId <= fEmcCrystals) { //calibrate as EMC
+ fADCchanelEmc = fCalibData->GetADCchannelEmc (module,column,row);
+ fADCpedestalEmc = fCalibData->GetADCpedestalEmc(module,column,row);
+ return fADCpedestalEmc + amp*fADCchanelEmc ;
+ }
+ else //calibrate as CPV, not implemented yet
+ return 0;
+ }
+ else{ //simulation
+ if(absId <= fEmcCrystals) //calibrate as EMC
+ return fADCpedestalEmc + amp*fADCchanelEmc ;
+ else //calibrate as CPV
+ return fADCpedestalCpv+ amp*fADCchanelCpv ;
+ }
+}
- if(!fIsInitialized) Init() ;
+//____________________________________________________________________________
+void AliPHOSClusterizerv1::Exec(Option_t *option)
+{
+ // Steering method to perform clusterization for events
+ // in the range from fFirstEvent to fLastEvent.
+ // This range is optionally set by SetEventRange().
+ // if fLastEvent=-1 (by default), then process events until the end.
if(strstr(option,"tim"))
- gBenchmark->Start("PHOSClusterizer");
+ gBenchmark->Start("PHOSClusterizer");
- Int_t nEvents = (Int_t) gAlice->TreeE()->GetEntries() ;
+ if(strstr(option,"print")) {
+ Print() ;
+ return ;
+ }
+
+ AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
+ if (fRawReader == 0)
+ gime->SetRawDigits(kFALSE);
+ else
+ gime->SetRawDigits(kTRUE);
+ GetCalibrationParameters() ;
- for(fEvent = 0;fEvent< nEvents; fEvent++){
- if(!ReadDigits()) //reads digits for event fEvent
- return;
- MakeClusters() ;
+ if (fLastEvent == -1)
+ fLastEvent = gime->MaxEvent() - 1 ;
+ else
+ fLastEvent = TMath::Min(fFirstEvent, gime->MaxEvent()); // one event at the time
+ Int_t nEvents = fLastEvent - fFirstEvent + 1;
+
+ Int_t ievent ;
+ for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) {
+ if (fRawReader == 0)
+ gime->Event(ievent ,"D"); // Read digits from simulated data
+ else
+ gime->Event(fRawReader,"W"); // Read digits from raw data
+
+ fNumberOfEmcClusters = fNumberOfCpvClusters = 0 ;
- if(fToUnfold) MakeUnfolding() ;
- WriteRecPoints() ;
- if(strstr(option,"deb"))
+ MakeClusters() ;
+
+ if(fToUnfold)
+ MakeUnfolding() ;
+
+ WriteRecPoints();
+
+ if(strstr(option,"deb"))
PrintRecPoints(option) ;
- }
+
+ //increment the total number of recpoints per run
+ fRecPointsInRun += gime->EmcRecPoints()->GetEntriesFast() ;
+ fRecPointsInRun += gime->CpvRecPoints()->GetEntriesFast() ;
+ }
+
+ if(fWrite) //do not unload in "on flight" mode
+ Unload();
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 ;
- }
-
-
+ AliInfo(Form("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
// The initial values for fitting procedure are set equal to the positions of local maxima.
// Cluster will be fitted as a superposition of nPar/3 electromagnetic showers
+
+ AliPHOSGetter * gime = AliPHOSGetter::Instance();
+ TClonesArray * digits = gime->Digits();
+
+
gMinuit->mncler(); // Reset Minuit's list of paramters
gMinuit->SetPrintLevel(-1) ; // No Printout
gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ;
TList * toMinuit = new TList();
toMinuit->AddAt(emcRP,0) ;
- toMinuit->AddAt(fDigits,1) ;
+ toMinuit->AddAt(digits,1) ;
gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
Int_t iDigit ;
+ 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 ;
- Float_t z ;
- fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
- fGeom->RelPosInModule(relid, x, z) ;
+ Float_t x = 0.;
+ Float_t z = 0.;
+ geom->AbsToRelNumbering(digit->GetId(), relid) ;
+ geom->RelPosInModule(relid, x, z) ;
Float_t energy = maxAtEnergy[iDigit] ;
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::Init()
+void AliPHOSClusterizerv1::GetCalibrationParameters()
{
- //Make all memory allocations which can not be done in default constructor.
- if(!fIsInitialized){
- if(fHeaderFileName.IsNull())
- fHeaderFileName = "galice.root" ;
-
-
- TFile * file = (TFile*) gROOT->GetFile(fHeaderFileName.Data() ) ;
-
- if(file == 0){
- file = new TFile(fHeaderFileName.Data(),"update") ;
- gAlice = (AliRun *) file->Get("gAlice") ;
- }
-
- AliPHOS * phos = (AliPHOS *) gAlice->GetDetector("PHOS") ;
- fGeom = AliPHOSGeometry::GetInstance(phos->GetGeometry()->GetName(),phos->GetGeometry()->GetTitle() );
-
- fDigits = new TClonesArray("AliPHOSDigit",10) ;
- fDigitizer = new AliPHOSDigitizer() ;
- fEmcRecPoints = new TObjArray(200) ;
- fCpvRecPoints = new TObjArray(200) ;
- fEmcRecPoints->SetOwner(); // This lets Clear() really detete rec.points in array
- fCpvRecPoints->SetOwner();
+ // Set calibration parameters:
+ // For raw data they are read from the calibration database,
+ // for simulated data they are taken from digitizer.
+ //
+ // It is a user responsilibity to open CDB and set
+ // AliPHOSCalibData object by the following operators:
+ //
+ // AliCDBLocal *loc = new AliCDBLocal("deCalibDB");
+ // AliPHOSCalibData* clb = (AliPHOSCalibData*)AliCDBStorage::Instance()
+ // ->Get(path_to_calibdata,run_number);
+ // AliPHOSGetter* gime = AliPHOSGetter::Instance("galice.root");
+ // gime->SetCalibData(clb);
+
+ AliPHOSGetter * gime = AliPHOSGetter::Instance();
+ if(gime->IsRawDigits()){
+ fCalibData = gime->CalibData();
+ }
+ else{
+ if ( !gime->Digitizer() )
+ gime->LoadDigitizer();
+ AliPHOSDigitizer * dig = gime->Digitizer();
+ fADCchanelEmc = dig->GetEMCchannel() ;
+ fADCpedestalEmc = dig->GetEMCpedestal();
- if(!gMinuit) gMinuit = new TMinuit(100) ;
+ fADCchanelCpv = dig->GetCPVchannel() ;
+ fADCpedestalCpv = dig->GetCPVpedestal() ;
+ }
+}
- // add Task to //root/Tasks folder
- TTask * roottasks = (TTask*)gROOT->GetRootFolder()->FindObject("Tasks") ;
- TTask * phostasks = (TTask*)(roottasks->GetListOfTasks()->FindObject("PHOS"));
- phostasks->Add(this) ;
+//____________________________________________________________________________
+void AliPHOSClusterizerv1::Init()
+{
+ // Make all memory allocations which can not be done in default constructor.
+ // Attach the Clusterizer task to the list of PHOS tasks
+
+ AliPHOSGetter* gime = AliPHOSGetter::Instance() ;
+ if(!gime)
+ gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data());
- fIsInitialized = kTRUE ;
+ AliPHOSGeometry * geom = gime->PHOSGeometry();
+
+ fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ;
+
+ if(!gMinuit)
+ gMinuit = new TMinuit(100);
+
+ if ( !gime->Clusterizer() ) {
+ gime->PostClusterizer(this);
}
}
+
+//____________________________________________________________________________
+void AliPHOSClusterizerv1::InitParameters()
+{
+
+ fNumberOfCpvClusters = 0 ;
+ fNumberOfEmcClusters = 0 ;
+
+ fCpvClusteringThreshold = 0.0;
+ fEmcClusteringThreshold = 0.2;
+
+ fEmcLocMaxCut = 0.03 ;
+ fCpvLocMaxCut = 0.03 ;
+
+ fEmcMinE = 0.01 ;
+ fCpvMinE = 0.0 ;
+
+ fW0 = 4.5 ;
+ fW0CPV = 4.0 ;
+
+ fEmcTimeGate = 1.e-8 ;
+
+ fToUnfold = kTRUE ;
+
+ fRecPointsInRun = 0 ;
+
+ fWrite = kTRUE ;
+
+ fCalibData = 0 ;
+
+ SetEventRange(0,-1) ;
+}
+
//____________________________________________________________________________
Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)const
{
// The order of d1 and d2 is important: first (d1) should be a digit already in a cluster
// which is compared to a digit (d2) not yet in a cluster
+ AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
+
Int_t rv = 0 ;
Int_t relid1[4] ;
- fGeom->AbsToRelNumbering(d1->GetId(), relid1) ;
+ geom->AbsToRelNumbering(d1->GetId(), relid1) ;
Int_t relid2[4] ;
- fGeom->AbsToRelNumbering(d2->GetId(), relid2) ;
+ 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 {
if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1))
- rv = 2; // Difference in row numbers is too large to look further
+ rv = 2; // Difference in row numbers is too large to look further
}
}
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] < fGeom->GetNumberOfPadsPhi()*fGeom->GetNumberOfPadsPhi() ) rv = 2 ;
-
return rv ;
}
+//____________________________________________________________________________
+void AliPHOSClusterizerv1::CleanDigits(TClonesArray * digits){
+ for(Int_t i=0; i<digits->GetEntriesFast(); i++){
+ AliPHOSDigit * digit = static_cast<AliPHOSDigit*>(digits->At(i)) ;
+ Float_t cut = IsInEmc(digit) ? fEmcMinE : fCpvMinE ;
+ if(Calibrate(digit->GetAmp(),digit->GetId()) < cut)
+ digits->RemoveAt(i) ;
+ }
+ digits->Compress() ;
+ for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) {
+ AliPHOSDigit *digit = static_cast<AliPHOSDigit*>( digits->At(i) ) ;
+ digit->SetIndexInList(i) ;
+ }
-
+}
//____________________________________________________________________________
Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const
{
// Tells if (true) or not (false) the digit is in a PHOS-EMC module
Bool_t rv = kFALSE ;
+ AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
- Int_t relid[4] ;
- fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
- if ( relid[1] == 0 ) rv = kTRUE;
+ if(digit->GetId() <= nEMC ) rv = kTRUE;
return rv ;
}
//____________________________________________________________________________
-Bool_t AliPHOSClusterizerv1::IsInPpsd(AliPHOSDigit * digit) const
+Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
{
- // Tells if (true) or not (false) the digit is in a PHOS-PPSD module
+ // Tells if (true) or not (false) the digit is in a PHOS-CPV module
Bool_t rv = kFALSE ;
+
+ AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
- Int_t relid[4] ;
- fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
- if ( relid[1] > 0 && relid[0] > fGeom->GetNCPVModules() ) rv = kTRUE;
+ if(digit->GetId() > nEMC ) rv = kTRUE;
return rv ;
}
//____________________________________________________________________________
-Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const
+void AliPHOSClusterizerv1::Unload()
{
- // Tells if (true) or not (false) the digit is in a PHOS-CPV module
-
- Bool_t rv = kFALSE ;
-
- Int_t relid[4] ;
- fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
-
- if ( relid[1] > 0 && relid[0] <= fGeom->GetNCPVModules() ) rv = kTRUE;
-
- return rv ;
+ AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
+ gime->PhosLoader()->UnloadDigits() ;
+ gime->PhosLoader()->UnloadRecPoints() ;
}
-//____________________________________________________________________________
-Bool_t AliPHOSClusterizerv1::ReadDigits()
- {
- // reads digitis with specified title from TreeD
-
- fNumberOfEmcClusters = 0 ;
- fNumberOfCpvClusters = 0 ;
- // Get Digits Tree header from file
- gAlice->GetEvent(fEvent) ;
- gAlice->SetEvent(fEvent) ;
-
- TTree * treeD = gAlice->TreeD() ;
-
- if(treeD==0){
- char treeName[20];
- sprintf(treeName,"TreeD%d",fEvent);
- cout << "Error in AliPHOSClusterizerv1 : no "<<treeName << endl ;
- cout << " Do nothing " << endl ;
- return kFALSE ;
- }
+//____________________________________________________________________________
+void AliPHOSClusterizerv1::WriteRecPoints()
+{
- TBranch * digitsBranch = 0;
- TBranch * digitizerBranch = 0;
+ // Creates new branches with given title
+ // fills and writes into TreeR.
- TObjArray * branches = treeD->GetListOfBranches() ;
- Int_t ibranch;
- Bool_t phosNotFound = kTRUE ;
- Bool_t digitizerNotFound = kTRUE ;
-
- for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
+ AliPHOSGetter * gime = AliPHOSGetter::Instance();
- if(phosNotFound){
- digitsBranch=(TBranch *) branches->At(ibranch) ;
- if( fDigitsBranchTitle.CompareTo(digitsBranch->GetTitle())==0 )
- if( strcmp(digitsBranch->GetName(),"PHOS") == 0)
- phosNotFound = kFALSE ;
- }
-
- if(digitizerNotFound){
- digitizerBranch = (TBranch *) branches->At(ibranch) ;
- if( fDigitsBranchTitle.CompareTo(digitizerBranch->GetTitle()) == 0)
- if( strcmp(digitizerBranch->GetName(),"AliPHOSDigitizer") == 0)
- digitizerNotFound = kFALSE ;
+ TObjArray * emcRecPoints = gime->EmcRecPoints() ;
+ TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
+ TClonesArray * digits = gime->Digits() ;
+
+
+ Int_t index ;
+ //Evaluate position, dispersion and other RecPoint properties..
+ Int_t nEmc = emcRecPoints->GetEntriesFast();
+ for(index = 0; index < nEmc; index++){
+ AliPHOSEmcRecPoint * rp = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) );
+ rp->Purify(fEmcMinE) ;
+ if(rp->GetMultiplicity()>0.) //If this RP is not empty
+ rp->EvalAll(fW0,digits) ;
+ else{
+ emcRecPoints->RemoveAt(index) ;
+ delete rp ;
}
-
}
-
- if(digitizerNotFound || phosNotFound){
- cout << "ERROR in AliPHOSClusterizerv1: " << endl ;
- cout << " Can't find Branch with digits or Digitizer "<< endl ; ;
- cout << " Do nothing" <<endl ;
- return kFALSE ;
+ emcRecPoints->Compress() ;
+ emcRecPoints->Sort() ;
+ // emcRecPoints->Expand(emcRecPoints->GetEntriesFast()) ;
+ for(index = 0; index < emcRecPoints->GetEntries(); index++){
+ dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) )->SetIndexInList(index) ;
}
- digitsBranch->SetAddress(&fDigits) ;
- digitizerBranch->SetAddress(&fDigitizer) ;
-
- digitsBranch->GetEntry(0) ;
- digitizerBranch->GetEntry(0) ;
-
- fPedestal = fDigitizer->GetPedestal() ;
- fSlope = fDigitizer->GetSlope() ;
- return kTRUE ;
-}
-
-//____________________________________________________________________________
-void AliPHOSClusterizerv1::WriteRecPoints()
-{
- // Checks, if PHOSEmcRP etc. branches with given title already exist,
- // if yes exits without writing,
- // else creates new branches with given title
- // fills and writes into TreeR.
-
- Int_t index ;
- //Evaluate poisition, dispersion and other RecPoint properties...
- for(index = 0; index < fEmcRecPoints->GetEntries(); index++)
- ((AliPHOSEmcRecPoint *)fEmcRecPoints->At(index))->EvalAll(fW0,fDigits) ;
-
- fEmcRecPoints->Sort() ;
-
- for(index = 0; index < fEmcRecPoints->GetEntries(); index++)
- ((AliPHOSEmcRecPoint *)fEmcRecPoints->At(index))->SetIndexInList(index) ;
-
- fEmcRecPoints->Expand(fEmcRecPoints->GetEntriesFast()) ;
-
//Now the same for CPV
- for(index = 0; index < fCpvRecPoints->GetEntries(); index++)
- ((AliPHOSRecPoint *)fCpvRecPoints->At(index))->EvalAll(fW0CPV,fDigits) ;
-
- fCpvRecPoints->Sort() ;
-
- for(index = 0; index < fCpvRecPoints->GetEntries(); index++)
- ((AliPHOSRecPoint *)fCpvRecPoints->At(index))->SetIndexInList(index) ;
-
- fCpvRecPoints->Expand(fCpvRecPoints->GetEntriesFast()) ;
-
- if(gAlice->TreeR()==0)
- gAlice->MakeTree("R") ;
+ for(index = 0; index < cpvRecPoints->GetEntries(); index++){
+ AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) );
+ rp->EvalAll(fW0CPV,digits) ;
+ }
+ cpvRecPoints->Sort() ;
-
- //Check, if branches already exist
- TBranch * emcBranch = 0;
- TBranch * cpvBranch = 0;
- TBranch * clusterizerBranch = 0;
+ for(index = 0; index < cpvRecPoints->GetEntries(); index++)
+ dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) )->SetIndexInList(index) ;
- TObjArray * branches = gAlice->TreeR()->GetListOfBranches() ;
- Int_t ibranch;
- Bool_t emcNotFound = kTRUE ;
- Bool_t cpvNotFound = kTRUE ;
- Bool_t clusterizerNotFound = kTRUE ;
+ cpvRecPoints->Expand(cpvRecPoints->GetEntriesFast()) ;
- for(ibranch = 0;ibranch <branches->GetEntries();ibranch++){
+ if(fWrite){ //We write TreeR
+ TTree * treeR = gime->TreeR();
- if(emcNotFound){
- emcBranch=(TBranch *) branches->At(ibranch) ;
- if( fRecPointsBranchTitle.CompareTo(emcBranch->GetTitle())==0 ){
- if( strcmp(emcBranch->GetName(),"PHOSEmcRP") == 0) {
- emcNotFound = kFALSE ;
- }
- }
- }
- if(cpvNotFound){
- cpvBranch=(TBranch *) branches->At(ibranch) ;
- if( fRecPointsBranchTitle.CompareTo(cpvBranch->GetTitle())==0 ){
- if( strcmp(cpvBranch->GetName(),"PHOSCpvRP") == 0) {
- cpvNotFound = kFALSE ;
- }
- }
- }
-
- if(clusterizerNotFound){
- clusterizerBranch = (TBranch *) branches->At(ibranch) ;
- if( fRecPointsBranchTitle.CompareTo(clusterizerBranch->GetTitle()) == 0){
- if( strcmp(clusterizerBranch->GetName(),"AliPHOSClusterizer") == 0) {
- clusterizerNotFound = kFALSE ;
- }
- }
- }
+ Int_t bufferSize = 32000 ;
+ Int_t splitlevel = 0 ;
- }
-
- if(!(clusterizerNotFound && emcNotFound && cpvNotFound)){
- cout << "AliPHOSClusterizer error" << endl;
- cout << " Branches PHOSEmcRP, PHOSCpvRP and AliPHOSClusterizer " << endl ;
- cout << " with title '" << fRecPointsBranchTitle.Data() <<"' already exist" << endl ;
- cout << " can not overwrite " << endl ;
- return ;
- }
-
- //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;
-
- //First EMC
- Int_t bufferSize = 32000 ;
- Int_t splitlevel = 0 ;
- emcBranch = gAlice->TreeR()->Branch("PHOSEmcRP","TObjArray",&fEmcRecPoints,bufferSize,splitlevel);
- emcBranch->SetTitle(fRecPointsBranchTitle.Data());
- 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
- cpvBranch = gAlice->TreeR()->Branch("PHOSCpvRP","TObjArray",&fCpvRecPoints,bufferSize,splitlevel);
- cpvBranch->SetTitle(fRecPointsBranchTitle.Data());
- if (filename) {
- cpvBranch->SetFile(filename);
- TIter next( cpvBranch->GetListOfBranches());
- TBranch * sb;
- while ((sb=(TBranch*)next())) {
- sb->SetFile(filename);
- }
- cwd->cd();
- }
+ //Now CPV branch
+ TBranch * cpvBranch = treeR->Branch("PHOSCpvRP","TObjArray",&cpvRecPoints,bufferSize,splitlevel);
+ cpvBranch->SetTitle(BranchName());
- //And Finally clusterizer branch
- AliPHOSClusterizerv1 * cl = this ;
- clusterizerBranch = gAlice->TreeR()->Branch("AliPHOSClusterizer","AliPHOSClusterizerv1",
- &cl,bufferSize,splitlevel);
- clusterizerBranch->SetTitle(fRecPointsBranchTitle.Data());
- if (filename) {
- clusterizerBranch->SetFile(filename);
- TIter next( clusterizerBranch->GetListOfBranches());
- TBranch * sb ;
- while ((sb=(TBranch*)next())) {
- sb->SetFile(filename);
- }
- cwd->cd();
+ emcBranch ->Fill() ;
+ cpvBranch ->Fill() ;
+
+ gime->WriteRecPoints("OVERWRITE");
+ gime->WriteClusterizer("OVERWRITE");
}
-
- emcBranch->Fill() ;
- cpvBranch->Fill() ;
- clusterizerBranch->Fill() ;
-
- gAlice->TreeR()->Write(0,kOverwrite) ;
-
}
//____________________________________________________________________________
{
// Steering method to construct the clusters stored in a list of Reconstructed Points
// A cluster is defined as a list of neighbour digits
- fEmcRecPoints->Clear() ;
- fCpvRecPoints->Clear() ;
+
+
+ AliPHOSGetter * gime = AliPHOSGetter::Instance();
+
+ TObjArray * emcRecPoints = gime->EmcRecPoints() ;
+ TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
+
+ emcRecPoints->Delete() ;
+ cpvRecPoints->Delete() ;
+ TClonesArray * digits = gime->Digits() ;
+
+ //Remove digits below threshold
+ CleanDigits(digits) ;
+
+ TClonesArray * digitsC = static_cast<TClonesArray*>( digits->Clone() ) ;
+
+
// Clusterization starts
- TClonesArray * digits = (TClonesArray*)fDigits->Clone() ;
- TIter nextdigit(digits) ;
+ TIter nextdigit(digitsC) ;
AliPHOSDigit * digit ;
Bool_t notremoved = kTRUE ;
- while ( (digit = (AliPHOSDigit *)nextdigit()) ) { // scan over the list of digits
+ while ( (digit = dynamic_cast<AliPHOSDigit *>( nextdigit()) ) ) { // scan over the list of digitsC
+
+
AliPHOSRecPoint * clu = 0 ;
TArrayI clusterdigitslist(1500) ;
Int_t index ;
- 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) ) {
- // start a new EMC RecPoint
- if(fNumberOfEmcClusters >= fEmcRecPoints->GetSize()) fEmcRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
- fEmcRecPoints->AddAt(new AliPHOSEmcRecPoint(), fNumberOfEmcClusters) ;
- clu = (AliPHOSEmcRecPoint *) fEmcRecPoints->At(fNumberOfEmcClusters) ;
- fNumberOfEmcClusters++ ;
- clu->AddDigit(*digit, Calibrate(digit->GetAmp())) ;
- clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
- iDigitInCluster++ ;
- digits->Remove(digit) ;
+ // start a new EMC RecPoint
+ if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
+ emcRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
+
+ emcRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
+ clu = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ;
+ fNumberOfEmcClusters++ ;
+ clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->GetId())) ;
+ clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
+ iDigitInCluster++ ;
+ digitsC->Remove(digit) ;
} else {
-
- // start a new PPSD/CPV cluster
- if(fNumberOfCpvClusters >= fCpvRecPoints->GetSize()) fCpvRecPoints->Expand(2*fNumberOfCpvClusters+1);
- if(IsInPpsd(digit))
- fCpvRecPoints->AddAt(new AliPHOSPpsdRecPoint(),fNumberOfCpvClusters) ;
- else
- fCpvRecPoints->AddAt(new AliPHOSCpvRecPoint(), fNumberOfCpvClusters) ;
- clu = (AliPHOSPpsdRecPoint *) fCpvRecPoints->At(fNumberOfCpvClusters) ;
- fNumberOfCpvClusters++ ;
-
- clu->AddDigit(*digit, Calibrate(digit->GetAmp()) ) ;
- clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
- iDigitInCluster++ ;
- digits->Remove(digit) ;
+
+ // start a new CPV cluster
+ if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
+ cpvRecPoints->Expand(2*fNumberOfCpvClusters+1);
+
+ cpvRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
+
+ clu = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;
+ fNumberOfCpvClusters++ ;
+ clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->GetId()) ) ;
+ clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
+ iDigitInCluster++ ;
+ digitsC->Remove(digit) ;
nextdigit.Reset() ;
-
- // Here we remove resting EMC digits, which cannot make cluster
-
+
+ // Here we remove remaining EMC digits, which cannot make a cluster
+
if( notremoved ) {
- while( ( digit = (AliPHOSDigit *)nextdigit() ) ) {
+ while( ( digit = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) {
if( IsInEmc(digit) )
- digits->Remove(digit) ;
+ digitsC->Remove(digit) ;
else
- break ;
- }
- notremoved = kFALSE ;
- }
-
+ break ;
+ }
+ notremoved = kFALSE ;
+ }
+
} // else
nextdigit.Reset() ;
AliPHOSDigit * digitN ;
index = 0 ;
while (index < iDigitInCluster){ // scan over digits already in cluster
- digit = (AliPHOSDigit*)fDigits->At(clusterdigitslist[index]) ;
- index++ ;
- while ( (digitN = (AliPHOSDigit *)nextdigit()) ) { // scan over the reduced list of digits
- Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
+ digit = dynamic_cast<AliPHOSDigit*>( digits->At(clusterdigitslist[index]) ) ;
+ index++ ;
+ while ( (digitN = dynamic_cast<AliPHOSDigit *>( nextdigit() ) ) ) { // scan over the reduced list of digits
+ Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!!
switch (ineb ) {
case 0 : // not a neighbour
- break ;
- case 1 : // are neighbours
- clu->AddDigit(*digitN, Calibrate( digitN->GetAmp() ) ) ;
- clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ;
- iDigitInCluster++ ;
- digits->Remove(digitN) ;
- break ;
+ break ;
+ case 1 : // are neighbours
+ clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), digitN->GetId() ) ) ;
+ clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ;
+ iDigitInCluster++ ;
+ digitsC->Remove(digitN) ;
+ break ;
case 2 : // too far from each other
- goto endofloop;
- } // switch
-
- } // while digitN
-
+ goto endofloop;
+ } // switch
+
+ } // while digitN
+
endofloop: ;
- nextdigit.Reset() ;
-
+ nextdigit.Reset() ;
+
} // loop over cluster
} // energy theshold
} // while digit
- delete digits ;
+ delete digitsC ;
}
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::Instance();
+ const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
+
+ TObjArray * emcRecPoints = gime->EmcRecPoints() ;
+ TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
+ TClonesArray * digits = gime->Digits() ;
+
// Unfold first EMC clusters
if(fNumberOfEmcClusters > 0){
- Int_t nModulesToUnfold = fGeom->GetNModules() ;
+ Int_t nModulesToUnfold = geom->GetNModules() ;
Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
Int_t index ;
for(index = 0 ; index < numberofNotUnfolded ; index++){
- AliPHOSEmcRecPoint * emcRecPoint = (AliPHOSEmcRecPoint *) fEmcRecPoints->At(index) ;
+ AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) ) ;
if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
- break ;
+ 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,fDigits) ;
+ 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) ;
- fEmcRecPoints->Remove(emcRecPoint);
- fEmcRecPoints->Compress() ;
- index-- ;
- fNumberOfEmcClusters -- ;
- numberofNotUnfolded-- ;
+ UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
+ emcRecPoints->Remove(emcRecPoint);
+ emcRecPoints->Compress() ;
+ index-- ;
+ fNumberOfEmcClusters -- ;
+ numberofNotUnfolded-- ;
+ }
+ else{
+ emcRecPoint->SetNExMax(1) ; //Only one local maximum
}
delete[] maxAt ;
// Unfold now CPV clusters
if(fNumberOfCpvClusters > 0){
- Int_t nModulesToUnfold = fGeom->GetNCPVModules() ;
+ Int_t nModulesToUnfold = geom->GetNModules() ;
Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
Int_t index ;
for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
- AliPHOSRecPoint * recPoint = (AliPHOSRecPoint *) fCpvRecPoints->At(index) ;
+ AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( cpvRecPoints->At(index) ) ;
if(recPoint->GetPHOSMod()> nModulesToUnfold)
- break ;
+ break ;
- AliPHOSEmcRecPoint * emcRecPoint = (AliPHOSEmcRecPoint*) recPoint ;
+ AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<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,fDigits) ;
+ Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,digits) ;
if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
- UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
- fCpvRecPoints->Remove(emcRecPoint);
- fCpvRecPoints->Compress() ;
- index-- ;
- numberofCpvNotUnfolded-- ;
- fNumberOfCpvClusters-- ;
+ UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
+ cpvRecPoints->Remove(emcRecPoint);
+ cpvRecPoints->Compress() ;
+ index-- ;
+ numberofCpvNotUnfolded-- ;
+ fNumberOfCpvClusters-- ;
}
delete[] maxAt ;
//____________________________________________________________________________
void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc,
- Int_t nMax,
- int * maxAt,
- Float_t * maxAtEnergy)
+ Int_t nMax,
+ AliPHOSDigit ** maxAt,
+ Float_t * maxAtEnergy)
{
// Performs the unfolding of a cluster with nMax overlapping showers
+ AliPHOSGetter * gime = AliPHOSGetter::Instance();
+
+ const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
+
+ const TClonesArray * digits = gime->Digits() ;
+ TObjArray * emcRecPoints = gime->EmcRecPoints() ;
+ TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
+
Int_t nPar = 3 * nMax ;
Float_t * fitparameters = new Float_t[nPar] ;
Bool_t rv = FindFit(iniEmc, maxAt, maxAtEnergy, nPar, fitparameters) ;
if( !rv ) {
// Fit failed, return and remove cluster
+ iniEmc->SetNExMax(-1) ;
delete[] fitparameters ;
return ;
}
Int_t nDigits = iniEmc->GetMultiplicity() ;
Float_t * efit = new Float_t[nDigits] ;
- Float_t xDigit,zDigit,distance ;
- Float_t xpar,zpar,epar ;
+ Float_t xDigit=0.,zDigit=0.,distance=0. ;
+ Float_t xpar=0.,zpar=0.,epar=0. ;
Int_t relid[4] ;
- AliPHOSDigit * digit ;
+ AliPHOSDigit * digit = 0 ;
Int_t * emcDigits = iniEmc->GetDigitsList() ;
Int_t iparam ;
Int_t iDigit ;
for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
- digit = (AliPHOSDigit*) fDigits->At(emcDigits[iDigit] ) ;
- fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
- fGeom->RelPosInModule(relid, xDigit, zDigit) ;
+ digit = dynamic_cast<AliPHOSDigit*>( digits->At(emcDigits[iDigit] ) ) ;
+ geom->AbsToRelNumbering(digit->GetId(), relid) ;
+ geom->RelPosInModule(relid, xDigit, zDigit) ;
efit[iDigit] = 0;
iparam = 0 ;
epar = fitparameters[iparam+2] ;
iparam += 3 ;
- AliPHOSEmcRecPoint * emcRP ;
+ AliPHOSEmcRecPoint * emcRP = 0 ;
if(iniEmc->IsEmc()){ //create new entries in fEmcRecPoints...
- if(fNumberOfEmcClusters >= fEmcRecPoints->GetSize())
- fEmcRecPoints->Expand(2*fNumberOfEmcClusters) ;
+ if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
+ emcRecPoints->Expand(2*fNumberOfEmcClusters) ;
- (*fEmcRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint() ;
- emcRP = (AliPHOSEmcRecPoint *) fEmcRecPoints->At(fNumberOfEmcClusters);
+ (*emcRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
+ emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ;
fNumberOfEmcClusters++ ;
+ emcRP->SetNExMax((Int_t)nPar/3) ;
}
else{//create new entries in fCpvRecPoints
- if(fNumberOfCpvClusters >= fCpvRecPoints->GetSize())
- fCpvRecPoints->Expand(2*fNumberOfCpvClusters) ;
+ if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
+ cpvRecPoints->Expand(2*fNumberOfCpvClusters) ;
- (*fCpvRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint() ;
- emcRP = (AliPHOSEmcRecPoint *) fCpvRecPoints->At(fNumberOfCpvClusters);
+ (*cpvRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
+ emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;
fNumberOfCpvClusters++ ;
}
Float_t eDigit ;
for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
- digit = (AliPHOSDigit*) fDigits->At( emcDigits[iDigit] ) ;
- fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
- fGeom->RelPosInModule(relid, xDigit, zDigit) ;
+ digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) ) ;
+ geom->AbsToRelNumbering(digit->GetId(), relid) ;
+ geom->RelPosInModule(relid, xDigit, zDigit) ;
distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ;
distance = TMath::Sqrt(distance) ;
ratio = epar * ShowerShape(distance) / efit[iDigit] ;
eDigit = emcEnergies[iDigit] * ratio ;
emcRP->AddDigit( *digit, eDigit ) ;
- }
+ }
}
delete[] fitparameters ;
// Calculates the Chi square for the cluster unfolding minimization
// Number of parameters, Gradient, Chi squared, parameters, what to do
+ TList * toMinuit = dynamic_cast<TList*>( gMinuit->GetObjectFit() ) ;
- TList * toMinuit = (TList*) gMinuit->GetObjectFit() ;
-
- AliPHOSEmcRecPoint * emcRP = (AliPHOSEmcRecPoint*) toMinuit->At(0) ;
- TClonesArray * digits = (TClonesArray*)toMinuit->At(1) ;
+ AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) ) ;
+ TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ;
- // AliPHOSEmcRecPoint * emcRP = (AliPHOSEmcRecPoint *) gMinuit->GetObjectFit() ; // EmcRecPoint to fit
+ // AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
Int_t * emcDigits = emcRP->GetDigitsList() ;
- Int_t nOfDigits = emcRP->GetDigitsMultiplicity() ;
+ Int_t nOdigits = emcRP->GetDigitsMultiplicity() ;
Float_t * emcEnergies = emcRP->GetEnergiesList() ;
-
- AliPHOSGeometry * geom = AliPHOSGeometry::GetInstance() ;
-
+
+ const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
fret = 0. ;
Int_t iparam ;
AliPHOSDigit * digit ;
Int_t iDigit ;
- for( iDigit = 0 ; iDigit < nOfDigits ; iDigit++) {
+ for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) {
- digit = (AliPHOSDigit*) digits->At( emcDigits[iDigit] ) ;
+ digit = dynamic_cast<AliPHOSDigit*>( digits->At( emcDigits[iDigit] ) );
Int_t relid[4] ;
Float_t xDigit ;
Int_t iParam = 0 ;
efit = 0 ;
while(iParam < nPar ){
- Double_t distance = (xDigit - x[iParam]) * (xDigit - x[iParam]) ;
- iParam++ ;
- distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ;
- distance = TMath::Sqrt( distance ) ;
- iParam++ ;
- efit += x[iParam] * ShowerShape(distance) ;
- iParam++ ;
+ Double_t distance = (xDigit - x[iParam]) * (xDigit - x[iParam]) ;
+ iParam++ ;
+ distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ;
+ distance = TMath::Sqrt( distance ) ;
+ iParam++ ;
+ efit += x[iParam] * ShowerShape(distance) ;
+ iParam++ ;
}
Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
iParam = 0 ;
while(iParam < nPar ){
- Double_t xpar = x[iParam] ;
- Double_t zpar = x[iParam+1] ;
- Double_t epar = x[iParam+2] ;
- Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
- Double_t shape = sum * ShowerShape(dr) ;
- Double_t r4 = dr*dr*dr*dr ;
- Double_t r295 = TMath::Power(dr,2.95) ;
- Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
- 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
-
- Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x
- iParam++ ;
- Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z
- iParam++ ;
- Grad[iParam] += shape ; // Derivative over energy
- iParam++ ;
+ Double_t xpar = x[iParam] ;
+ Double_t zpar = x[iParam+1] ;
+ Double_t epar = x[iParam+2] ;
+ Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
+ Double_t shape = sum * ShowerShape(dr) ;
+ Double_t r4 = dr*dr*dr*dr ;
+ Double_t r295 = TMath::Power(dr,2.95) ;
+ Double_t deriv =-4. * dr*dr * ( 2.32 / ( (2.32 + 0.26 * r4) * (2.32 + 0.26 * r4) ) +
+ 0.0316 * (1. + 0.0171 * r295) / ( ( 1. + 0.0652 * r295) * (1. + 0.0652 * r295) ) ) ;
+
+ Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x
+ iParam++ ;
+ Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z
+ iParam++ ;
+ Grad[iParam] += shape ; // Derivative over energy
+ iParam++ ;
}
}
efit = 0;
}
//____________________________________________________________________________
-void AliPHOSClusterizerv1::Print(Option_t * option)const
+void AliPHOSClusterizerv1::Print(const Option_t *)const
{
// Print clusterizer parameters
- if(fIsInitialized){
-
+ TString message ;
+ TString taskName(GetName()) ;
+ taskName.ReplaceAll(Version(), "") ;
+
+ if( strcmp(GetName(), "") !=0 ) {
// Print parameters
-
- cout << "---------------"<< GetName() << " " << 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 " ;
+
+ AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", 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
- cout << "AliPHOSClusterizerv1: " << endl ;
- cout << " Found "<< fEmcRecPoints->GetEntriesFast() << " EMC Rec Points and "
- << fCpvRecPoints->GetEntriesFast() << " CPV RecPoints" << endl ;
+ AliPHOSGetter * gime = AliPHOSGetter::Instance();
+
+ TObjArray * emcRecPoints = gime->EmcRecPoints() ;
+ TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
+ AliInfo(Form("\nevent %d \n Found %d EMC RecPoints and %d CPV RecPoints",
+ gAlice->GetEvNumber(),
+ emcRecPoints->GetEntriesFast(),
+ cpvRecPoints->GetEntriesFast() )) ;
+
+ fRecPointsInRun += emcRecPoints->GetEntriesFast() ;
+ fRecPointsInRun += cpvRecPoints->GetEntriesFast() ;
+
+
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;
-
+ printf("\n EMC clusters \n") ;
+ printf("Index Ene(MeV) Multi Module X Y Z Lambdas_1 Lambda_2 # of prim Primaries list\n") ;
Int_t index ;
- for (index = 0 ; index < fEmcRecPoints->GetEntries() ; index++) {
- AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEmcRecPoints->At(index) ;
-
- cout << setw(6) << rp->GetIndexInList() << " ";
- cout << setw(6) << rp->GetEnergy() << " ";
- cout << setw(6) << rp->GetMultiplicity()<< " ";
- cout << setw(6) << rp->GetPHOSMod() << " ";
-
+ for (index = 0 ; index < emcRecPoints->GetEntries() ; index++) {
+ AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )emcRecPoints->At(index) ;
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 PCV/PPSD recPoints
- cout << "EMC clusters " << endl ;
- cout << " Index "
- << " Multi "
- << " Module "
- << " Layer "
- << " X "
- << " Y "
- << " Z "
- << " # of prim "
- << " Primaries list " << endl;
-
- for (index = 0 ; index < fEmcRecPoints->GetEntries() ; index++) {
- AliPHOSRecPoint * rp = (AliPHOSRecPoint * )fCpvRecPoints->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 ";
+ 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] ) ;
}
- else
- cout <<" CPV ";
+ printf("\n") ;
+ }
+
+ //Now plot CPV recPoints
+ printf("\n CPV clusters \n") ;
+ printf("Index Ene(MeV) Module X Y Z \n") ;
+ for (index = 0 ; index < cpvRecPoints->GetEntries() ; index++) {
+ AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )cpvRecPoints->At(index) ;
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;
+ printf("\n%6d %8.2f %2d %4.1f %4.1f %4.1f \n",
+ rp->GetIndexInList(), rp->GetEnergy(), rp->GetPHOSMod(),
+ locpos.X(), locpos.Y(), locpos.Z()) ;
}
-
-
- cout << "-------------------------------------------------"<<endl ;
}
}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff