/* $Id$ */
+/* History of cvs commits:
+ *
+ * $Log$
+ * Revision 1.93 2006/08/01 12:20:17 cvetan
+ * 1. Adding a possibility to read and reconstruct an old rcu formatted raw data. This is controlled by an option of AliReconstruction and AliPHOSReconstructor. 2. In case of raw data processing (without galice.root) create the default AliPHOSGeometry object. Most likely this should be moved to the CDB
+ *
+ * Revision 1.92 2006/04/29 20:26:46 hristov
+ * Separate EMC and CPV calibration (Yu.Kharlov)
+ *
+ * Revision 1.91 2006/04/22 10:30:17 hristov
+ * Add fEnergy to AliPHOSDigit and operate with EMC amplitude in energy units (Yu.Kharlov)
+ *
+ * Revision 1.90 2006/04/11 15:22:59 hristov
+ * run number in query set to -1: forces AliCDBManager to use its run number (A.Colla)
+ *
+ * Revision 1.89 2006/03/13 14:05:42 kharlov
+ * Calibration objects for EMC and CPV
+ *
+ * Revision 1.88 2006/01/11 08:54:52 hristov
+ * Additional protection in case no calibration entry was found
+ *
+ * Revision 1.87 2005/11/22 08:46:43 kharlov
+ * Updated to new CDB (Boris Polichtchouk)
+ *
+ * Revision 1.86 2005/11/14 21:52:43 hristov
+ * Coding conventions
+ *
+ * Revision 1.85 2005/09/27 16:08:08 hristov
+ * New version of CDB storage framework (A.Colla)
+ *
+ * Revision 1.84 2005/09/21 10:02:47 kharlov
+ * Reading calibration from CDB (Boris Polichtchouk)
+ *
+ * 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
// --- Standard library ---
// --- AliRoot header files ---
+#include "AliLog.h"
#include "AliPHOSGetter.h"
#include "AliPHOSGeometry.h"
#include "AliPHOSClusterizerv1.h"
#include "AliPHOSCpvRecPoint.h"
#include "AliPHOSDigit.h"
#include "AliPHOSDigitizer.h"
+#include "AliPHOSCalibrationDB.h"
+#include "AliCDBManager.h"
+#include "AliCDBStorage.h"
+#include "AliCDBEntry.h"
ClassImp(AliPHOSClusterizerv1)
{
// dtor
- // delete fPedestals ;
-// delete fGains ;
}
-
-
//____________________________________________________________________________
const TString AliPHOSClusterizerv1::BranchName() const
{
}
//____________________________________________________________________________
-Float_t AliPHOSClusterizerv1::Calibrate(Int_t amp, Int_t absId) const
+Float_t AliPHOSClusterizerv1::CalibrateEMC(Float_t amp, Int_t absId)
{
- //To be replaced later by the method, reading individual parameters from the database
- // if(fCalibrationDB)
- // return fCalibrationDB->Calibrate(amp,absId) ;
- // else{ //simulation
- if(absId <= fEmcCrystals) //calibrate as EMC
- return fADCpedestalEmc + amp*fADCchanelEmc ;
- else //calibrate as CPV
- return fADCpedestalCpv+ amp*fADCchanelCpv ;
- // }
+ // Convert EMC measured amplitude into real 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 column = relId[3];
+ Int_t row = relId[2];
+ if(absId <= fEmcCrystals) { // this is EMC
+ fADCchanelEmc = fCalibData->GetADCchannelEmc (module,column,row);
+ return amp*fADCchanelEmc ;
+ }
+ }
+ else{ //simulation
+ if(absId <= fEmcCrystals) // this is EMC
+ return fADCpedestalEmc + amp*fADCchanelEmc ;
+ }
+ return 0;
+}
+
+//____________________________________________________________________________
+Float_t AliPHOSClusterizerv1::CalibrateCPV(Int_t amp, Int_t absId)
+{
+ // Convert digitized CPV amplitude into charge.
+ // 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 column = relId[3];
+ Int_t row = relId[2];
+ if(absId > fEmcCrystals) { // this is CPV
+ fADCchanelCpv = fCalibData->GetADCchannelCpv (module,column,row);
+ fADCpedestalCpv = fCalibData->GetADCpedestalCpv(module,column,row);
+ return fADCpedestalCpv + amp*fADCchanelCpv ;
+ }
+ }
+ else{ //simulation
+ if(absId > fEmcCrystals) // this is CPV
+ return fADCpedestalCpv+ amp*fADCchanelCpv ;
+ }
+ return 0;
}
//____________________________________________________________________________
return ;
}
- AliPHOSGetter * gime = AliPHOSGetter::Instance(GetTitle()) ;
+ GetCalibrationParameters() ;
+
+ AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
+ if (fRawReader == 0)
+ gime->SetRawDigits(kFALSE);
+ else
+ gime->SetRawDigits(kTRUE);
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++) {
- gime->Event(ievent, "D");
-
- GetCalibrationParameters() ;
-
+ if (fRawReader == 0)
+ gime->Event(ievent ,"D"); // Read digits from simulated data
+ else {
+ gime->Event(fRawReader,"W",fIsOldRCUFormat); // Read digits from raw data
+ }
fNumberOfEmcClusters = fNumberOfCpvClusters = 0 ;
MakeClusters() ;
//increment the total number of recpoints per run
fRecPointsInRun += gime->EmcRecPoints()->GetEntriesFast() ;
fRecPointsInRun += gime->CpvRecPoints()->GetEntriesFast() ;
- }
+ }
- Unload();
+ if(fWrite) //do not unload in "on flight" mode
+ Unload();
if(strstr(option,"tim")){
gBenchmark->Stop("PHOSClusterizer");
- Info("Exec", " took %f seconds for Clusterizing %f seconds per event \n",
+ AliInfo(Form("took %f seconds for Clusterizing %f seconds per event \n",
gBenchmark->GetCpuTime("PHOSClusterizer"),
- gBenchmark->GetCpuTime("PHOSClusterizer")/nEvents ) ;
+ gBenchmark->GetCpuTime("PHOSClusterizer")/nEvents )) ;
}
}
//____________________________________________________________________________
void AliPHOSClusterizerv1::GetCalibrationParameters()
{
+ // Set calibration parameters:
+ // if calibration database exists, they are read from database,
+ // otherwise, they are taken from digitizer.
+ //
+ // It is a user responsilibity to open CDB before reconstruction, for example:
+ // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
AliPHOSGetter * gime = AliPHOSGetter::Instance();
-
- if ( !gime->Digitizer() )
- gime->LoadDigitizer();
- AliPHOSDigitizer * dig = gime->Digitizer();
- fADCchanelEmc = dig->GetEMCchannel() ;
- fADCpedestalEmc = dig->GetEMCpedestal();
+ // fCalibData = new AliPHOSCalibData(gAlice->GetRunNumber()); //original
+
+ fCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
- fADCchanelCpv = dig->GetCPVchannel() ;
- fADCpedestalCpv = dig->GetCPVpedestal() ;
-
- // fCalibrationDB = gime->CalibrationDB();
+ if(!fCalibData)
+ {
+ if ( !gime->Digitizer() )
+ gime->LoadDigitizer();
+ AliPHOSDigitizer * dig = gime->Digitizer();
+ fADCchanelEmc = dig->GetEMCchannel() ;
+ fADCpedestalEmc = dig->GetEMCpedestal();
+
+ fADCchanelCpv = dig->GetCPVchannel() ;
+ fADCpedestalCpv = dig->GetCPVpedestal() ;
+ }
}
//____________________________________________________________________________
{
// 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(GetTitle(), fEventFolderName.Data());
+
+ AliPHOSGetter* gime = AliPHOSGetter::Instance() ;
+ if(!gime)
+ gime = AliPHOSGetter::Instance(GetTitle(), fEventFolderName.Data());
AliPHOSGeometry * geom = gime->PHOSGeometry();
-
+ if (!geom) {
+ AliError("Could not find PHOS geometry! Loading the default one !");
+ geom = AliPHOSGeometry::GetInstance("IHEP","");
+ }
+
fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ;
if(!gMinuit)
gMinuit = new TMinuit(100);
- if ( !gime->Clusterizer() )
+ if ( !gime->Clusterizer() ) {
gime->PostClusterizer(this);
+ }
}
//____________________________________________________________________________
fEmcLocMaxCut = 0.03 ;
fCpvLocMaxCut = 0.03 ;
-
+
+ fEmcMinE = 0.01 ;
+ fCpvMinE = 0.0 ;
+
fW0 = 4.5 ;
fW0CPV = 4.0 ;
fRecPointsInRun = 0 ;
+ fWrite = kTRUE ;
+
+ fCalibData = 0 ;
+
SetEventRange(0,-1) ;
+
+ fIsOldRCUFormat = kFALSE;
}
//____________________________________________________________________________
return rv ;
}
+//____________________________________________________________________________
+void AliPHOSClusterizerv1::CleanDigits(TClonesArray * digits)
+{
+ // Remove digits with amplitudes below threshold
+ for(Int_t i=0; i<digits->GetEntriesFast(); i++){
+ AliPHOSDigit * digit = static_cast<AliPHOSDigit*>(digits->At(i)) ;
+ if ( (IsInEmc(digit) && CalibrateEMC(digit->GetEnergy(),digit->GetId()) < fEmcMinE) ||
+ (IsInCpv(digit) && CalibrateCPV(digit->GetAmp() ,digit->GetId()) < fCpvMinE) )
+ 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
{
TObjArray * emcRecPoints = gime->EmcRecPoints() ;
TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
TClonesArray * digits = gime->Digits() ;
-
- TTree * treeR = gime->TreeR();
-
+
+
Int_t index ;
- //Evaluate position, dispersion and other RecPoint properties...
- for(index = 0; index < emcRecPoints->GetEntries(); index++)
- dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) )->EvalAll(fW0,digits) ;
-
+ //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 ;
+ }
+ }
+ emcRecPoints->Compress() ;
emcRecPoints->Sort() ;
- for(index = 0; index < emcRecPoints->GetEntries(); index++)
+ // emcRecPoints->Expand(emcRecPoints->GetEntriesFast()) ;
+ for(index = 0; index < emcRecPoints->GetEntries(); index++){
dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) )->SetIndexInList(index) ;
-
- emcRecPoints->Expand(emcRecPoints->GetEntriesFast()) ;
+ }
//Now the same for CPV
- for(index = 0; index < cpvRecPoints->GetEntries(); index++)
- dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) )->EvalAll(fW0CPV,digits) ;
-
+ for(index = 0; index < cpvRecPoints->GetEntries(); index++){
+ AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) );
+ rp->EvalAll(fW0CPV,digits) ;
+ }
cpvRecPoints->Sort() ;
for(index = 0; index < cpvRecPoints->GetEntries(); index++)
cpvRecPoints->Expand(cpvRecPoints->GetEntriesFast()) ;
- Int_t bufferSize = 32000 ;
- Int_t splitlevel = 0 ;
-
- //First EMC
- TBranch * emcBranch = treeR->Branch("PHOSEmcRP","TObjArray",&emcRecPoints,bufferSize,splitlevel);
- emcBranch->SetTitle(BranchName());
+ if(fWrite){ //We write TreeR
+ TTree * treeR = gime->TreeR();
- //Now CPV branch
- TBranch * cpvBranch = treeR->Branch("PHOSCpvRP","TObjArray",&cpvRecPoints,bufferSize,splitlevel);
- cpvBranch->SetTitle(BranchName());
+ Int_t bufferSize = 32000 ;
+ Int_t splitlevel = 0 ;
- emcBranch ->Fill() ;
- cpvBranch ->Fill() ;
-
- gime->WriteRecPoints("OVERWRITE");
- gime->WriteClusterizer("OVERWRITE");
+ //First EMC
+ TBranch * emcBranch = treeR->Branch("PHOSEmcRP","TObjArray",&emcRecPoints,bufferSize,splitlevel);
+ emcBranch->SetTitle(BranchName());
+
+ //Now CPV branch
+ TBranch * cpvBranch = treeR->Branch("PHOSCpvRP","TObjArray",&cpvRecPoints,bufferSize,splitlevel);
+ cpvBranch->SetTitle(BranchName());
+
+ emcBranch ->Fill() ;
+ cpvBranch ->Fill() ;
+
+ gime->WriteRecPoints("OVERWRITE");
+ gime->WriteClusterizer("OVERWRITE");
+ }
}
//____________________________________________________________________________
TClonesArray * digits = gime->Digits() ;
+ //Remove digits below threshold
+ CleanDigits(digits) ;
+
TClonesArray * digitsC = static_cast<TClonesArray*>( digits->Clone() ) ;
-
// Clusterization starts
TArrayI clusterdigitslist(1500) ;
Int_t index ;
- if (( IsInEmc (digit) && Calibrate(digit->GetAmp(),digit->GetId()) > fEmcClusteringThreshold ) ||
- ( IsInCpv (digit) && Calibrate(digit->GetAmp(),digit->GetId()) > fCpvClusteringThreshold ) ) {
+ if (( IsInEmc (digit) &&
+ CalibrateEMC(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) ||
+ ( IsInCpv (digit) &&
+ CalibrateCPV(digit->GetAmp() ,digit->GetId()) > fCpvClusteringThreshold ) ) {
Int_t iDigitInCluster = 0 ;
if ( IsInEmc(digit) ) {
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++ ;
+ fNumberOfEmcClusters++ ;
+ clu->AddDigit(*digit, CalibrateEMC(digit->GetEnergy(),digit->GetId())) ;
+ clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
+ iDigitInCluster++ ;
digitsC->Remove(digit) ;
} else {
clu = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;
fNumberOfCpvClusters++ ;
- clu->AddDigit(*digit, Calibrate(digit->GetAmp(),digit->GetId()) ) ;
+ clu->AddDigit(*digit, CalibrateCPV(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(), digitN->GetId() ) ) ;
+ if (IsInEmc (digitN))
+ clu->AddDigit(*digitN, CalibrateEMC( digitN->GetEnergy(), digitN->GetId() ) );
+ else
+ clu->AddDigit(*digitN, CalibrateCPV( digitN->GetAmp() , digitN->GetId() ) );
+
clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ;
iDigitInCluster++ ;
digitsC->Remove(digitN) ;
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
+ iniEmc->SetNExMax(-1) ;
delete[] fitparameters ;
return ;
}
(*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 >= cpvRecPoints->GetSize())
}
//____________________________________________________________________________
-void AliPHOSClusterizerv1::Print()const
+void AliPHOSClusterizerv1::Print(const Option_t *)const
{
// Print clusterizer parameters
else
message = " AliPHOSClusterizerv1 not initialized " ;
- Info("Print", message.Data(),
+ AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(),
taskName.Data(),
GetTitle(),
taskName.Data(),
fW0,
fCpvClusteringThreshold,
fCpvLocMaxCut,
- fW0CPV ) ;
+ fW0CPV )) ;
}
TObjArray * emcRecPoints = gime->EmcRecPoints() ;
TObjArray * cpvRecPoints = gime->CpvRecPoints() ;
- printf("\nevent %d \n", gAlice->GetEvNumber()) ;
- printf(" Found %d EMC RecPoints and %d CPV RecPoints \n",
- emcRecPoints->GetEntriesFast(),cpvRecPoints->GetEntriesFast()) ;
+ 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() ;