/* History of cvs commits:
*
* $Log$
+ * Revision 1.112 2007/08/22 09:20:50 hristov
+ * Updated QA classes (Yves)
+ *
+ * Revision 1.111 2007/08/08 12:11:28 kharlov
+ * Protection against uninitialized fQADM
+ *
+ * Revision 1.110 2007/08/07 14:16:00 kharlov
+ * Quality assurance added (Yves Schutz)
+ *
+ * Revision 1.109 2007/07/24 17:20:35 policheh
+ * Usage of RecoParam objects instead of hardcoded parameters in reconstruction.
+ * (See $ALICE_ROOT/PHOS/macros/BeamTest2006/RawReconstruction.C).
+ *
+ * Revision 1.108 2007/06/18 07:00:51 kharlov
+ * Bug fix for attempt to use AliPHOSEmcRecPoint after its deletion
+ *
+ * Revision 1.107 2007/05/25 14:12:26 policheh
+ * Local to tracking CS transformation added for CPV rec. points
+ *
+ * Revision 1.106 2007/05/24 13:01:22 policheh
+ * Local to tracking CS transformation invoked for each EMC rec.point
+ *
+ * Revision 1.105 2007/05/02 13:41:22 kharlov
+ * Mode protection against absence of calib.data from AliPHOSCalibData to AliPHOSClusterizerv1::GetCalibrationParameters()
+ *
+ * Revision 1.104 2007/04/27 16:55:53 kharlov
+ * Calibration stops if PHOS CDB objects do not exist
+ *
+ * Revision 1.103 2007/04/11 11:55:45 policheh
+ * SetDistancesToBadChannels() added.
+ *
+ * Revision 1.102 2007/03/28 19:18:15 kharlov
+ * RecPoints recalculation in TSM removed
+ *
+ * Revision 1.101 2007/03/06 06:51:27 kharlov
+ * Calculation of cluster properties dep. on vertex posponed to TrackSegmentMaker
+ *
+ * Revision 1.100 2007/01/10 11:07:26 kharlov
+ * Raw digits writing to file (B.Polichtchouk)
+ *
+ * Revision 1.99 2006/11/07 16:49:51 kharlov
+ * Corrections for next event switch in case of raw data (B.Polichtchouk)
+ *
+ * Revision 1.98 2006/10/27 17:14:27 kharlov
+ * Introduce AliDebug and AliLog (B.Polichtchouk)
+ *
* Revision 1.97 2006/08/29 11:41:19 kharlov
* Missing implementation of ctors and = operator are added
*
// 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", "recpointsname", "digitsname")
-// Warning in <TDatabasePDG::TDatabasePDG>: object already instantiated
-// // 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 [0] AliPHOSClusterizerv1 * cl = new AliPHOSClusterizerv1(<pointer_to_phos_geometry_onject>)
+// root [1] cl->Digits2Clusters(digitsTree,clusterTree)
+// //finds RecPoints in the current event
// root [2] cl->SetDigitsBranch("digits2")
// //sets another title for Digitis (input) branch
// root [3] cl->SetRecPointsBranch("recp2")
// //sets another title four output branches
// root [4] cl->SetEmcLocalMaxCut(0.03)
// //set clusterization parameters
-// root [5] cl->ExecuteTask("deb all time")
-// //once more finds RecPoints options are
-// // deb - print number of found rec points
-// // deb all - print number of found RecPoints and some their characteristics
-// // time - print benchmarking results
// --- ROOT system ---
#include "TMinuit.h"
#include "TTree.h"
#include "TBenchmark.h"
+#include "TClonesArray.h"
// --- Standard library ---
// --- AliRoot header files ---
#include "AliLog.h"
-#include "AliPHOSGetter.h"
+#include "AliConfig.h"
#include "AliPHOSGeometry.h"
#include "AliPHOSClusterizerv1.h"
#include "AliPHOSEmcRecPoint.h"
#include "AliCDBManager.h"
#include "AliCDBStorage.h"
#include "AliCDBEntry.h"
+#include "AliPHOSRecoParam.h"
+#include "AliPHOSQualAssDataMaker.h"
+#include "AliPHOSCalibData.h"
+#include "AliPHOSReconstructor.h"
ClassImp(AliPHOSClusterizerv1)
fADCchanelCpv(0), fADCpedestalCpv(0), fEmcClusteringThreshold(0),
fCpvClusteringThreshold(0), fEmcMinE(0), fCpvMinE(0),
fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
- fW0CPV(0), fRecPointsInRun(0), fEmcTimeGate(0),
+ fW0CPV(0), fEmcTimeGate(0),
fIsOldRCUFormat(0)
{
// default ctor (to be used mainly by Streamer)
}
//____________________________________________________________________________
-AliPHOSClusterizerv1::AliPHOSClusterizerv1(const TString alirunFileName, const TString eventFolderName) :
- AliPHOSClusterizer(alirunFileName, eventFolderName),
+AliPHOSClusterizerv1::AliPHOSClusterizerv1(AliPHOSGeometry *geom) :
+ AliPHOSClusterizer(geom),
fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
fWrite(0), fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
fCalibData(0), fADCchanelEmc(0), fADCpedestalEmc(0),
fADCchanelCpv(0), fADCpedestalCpv(0), fEmcClusteringThreshold(0),
fCpvClusteringThreshold(0), fEmcMinE(0), fCpvMinE(0),
fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
- fW0CPV(0), fRecPointsInRun(0), fEmcTimeGate(0),
+ fW0CPV(0), fEmcTimeGate(0),
fIsOldRCUFormat(0)
{
// ctor with the indication of the file where header Tree and digits Tree are stored
fDefaultInit = kFALSE ;
}
-//____________________________________________________________________________
-AliPHOSClusterizerv1::AliPHOSClusterizerv1(const AliPHOSClusterizerv1 & obj) :
- AliPHOSClusterizer(obj),
- fDefaultInit(0), fEmcCrystals(0), fToUnfold(0),
- fWrite(0), fNumberOfEmcClusters(0), fNumberOfCpvClusters(0),
- fCalibData(0), fADCchanelEmc(0), fADCpedestalEmc(0),
- fADCchanelCpv(0), fADCpedestalCpv(0), fEmcClusteringThreshold(0),
- fCpvClusteringThreshold(0), fEmcMinE(0), fCpvMinE(0),
- fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0),
- fW0CPV(0), fRecPointsInRun(0), fEmcTimeGate(0),
- fIsOldRCUFormat(0)
-{
- // Copy constructor
-}
//____________________________________________________________________________
AliPHOSClusterizerv1::~AliPHOSClusterizerv1()
{
// dtor
}
-//____________________________________________________________________________
-const TString AliPHOSClusterizerv1::BranchName() const
-{
- return GetName();
-}
-
+
//____________________________________________________________________________
Float_t AliPHOSClusterizerv1::CalibrateEMC(Float_t amp, Int_t absId)
{
if(fCalibData){
Int_t relId[4];
- AliPHOSGetter *gime = AliPHOSGetter::Instance();
- gime->PHOSGeometry()->AbsToRelNumbering(absId,relId) ;
+ fGeom->AbsToRelNumbering(absId,relId) ;
Int_t module = relId[0];
Int_t column = relId[3];
Int_t row = relId[2];
if(fCalibData){
Int_t relId[4];
- AliPHOSGetter *gime = AliPHOSGetter::Instance();
- gime->PHOSGeometry()->AbsToRelNumbering(absId,relId) ;
+ fGeom->AbsToRelNumbering(absId,relId) ;
Int_t module = relId[0];
Int_t column = relId[3];
Int_t row = relId[2];
}
//____________________________________________________________________________
-void AliPHOSClusterizerv1::Exec(Option_t *option)
+void AliPHOSClusterizerv1::Digits2Clusters(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.
+ // Steering method to perform clusterization for one event
+ // The input is the tree with digits.
+ // The output is the tree with clusters.
if(strstr(option,"tim"))
gBenchmark->Start("PHOSClusterizer");
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++) {
- 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() ;
+ MakeClusters() ;
- AliDebug(2,Form(" ---- Printing clusters (%d) of event %d.\n",
- gime->EmcRecPoints()->GetEntries(),ievent));
- if(AliLog::GetGlobalDebugLevel()>1)
- gime->EmcRecPoints()->Print();
+ AliDebug(2,Form(" ---- Printing clusters (%d)\n",
+ fEMCRecPoints->GetEntries()));
+ if(AliLog::GetGlobalDebugLevel()>1)
+ fEMCRecPoints->Print();
+
+ if(fToUnfold)
+ MakeUnfolding();
+
+ //makes the quality assurance data
+ if (GetQualAssDataMaker()) {
+ GetQualAssDataMaker()->SetData(fEMCRecPoints) ;
+ GetQualAssDataMaker()->Exec(AliQualAss::kRECPOINTS) ;
+ GetQualAssDataMaker()->SetData(fCPVRecPoints) ;
+ GetQualAssDataMaker()->Exec(AliQualAss::kRECPOINTS) ;
+ }
- if(fToUnfold)
- MakeUnfolding() ;
+ WriteRecPoints();
- WriteRecPoints();
+ if(strstr(option,"deb"))
+ PrintRecPoints(option) ;
- if(strstr(option,"deb"))
- PrintRecPoints(option) ;
+ // PLEASE FIX BY MOVING IT TO ALIRECONSTRUCTION !!!
+ //Write the quality assurance data only after the last event
+ // if (GetQualAssDataMaker() && fEventCounter == gime->MaxEvent())
+ // GetQualAssDataMaker()->Finish(AliQualAss::kRECPOINTS) ;
- //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");
- AliInfo(Form("took %f seconds for Clusterizing %f seconds per event \n",
- gBenchmark->GetCpuTime("PHOSClusterizer"),
- gBenchmark->GetCpuTime("PHOSClusterizer")/nEvents )) ;
+ AliInfo(Form("took %f seconds for Clusterizing\n",
+ gBenchmark->GetCpuTime("PHOSClusterizer")));
}
}
// 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(digits,1) ;
+ toMinuit->AddAt(fDigitsArr,1) ;
+ toMinuit->AddAt(fGeom,2) ;
gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare
Int_t iDigit ;
- const AliPHOSGeometry * geom = gime->PHOSGeometry() ;
-
for(iDigit = 0; iDigit < nDigits; iDigit++){
digit = maxAt[iDigit];
Int_t relid[4] ;
Float_t x = 0.;
Float_t z = 0.;
- geom->AbsToRelNumbering(digit->GetId(), relid) ;
- geom->RelPosInModule(relid, x, z) ;
+ fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ fGeom->RelPosInModule(relid, x, z) ;
Float_t energy = maxAtEnergy[iDigit] ;
{
// Set calibration parameters:
// if calibration database exists, they are read from database,
- // otherwise, they are taken from digitizer.
+ // otherwise, reconstruction stops in the constructor of AliPHOSCalibData
//
// It is a user responsilibity to open CDB before reconstruction, for example:
// AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
- AliPHOSGetter * gime = AliPHOSGetter::Instance();
- // fCalibData = new AliPHOSCalibData(gAlice->GetRunNumber()); //original
-
fCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number
-
- if(!fCalibData)
- {
- if ( !gime->Digitizer() )
- gime->LoadDigitizer();
- AliPHOSDigitizer * dig = gime->Digitizer();
- fADCchanelEmc = dig->GetEMCchannel() ;
- fADCpedestalEmc = dig->GetEMCpedestal();
-
- fADCchanelCpv = dig->GetCPVchannel() ;
- fADCpedestalCpv = dig->GetCPVpedestal() ;
- }
+ if (fCalibData->GetCalibDataEmc() == 0)
+ AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n");
+ if (fCalibData->GetCalibDataCpv() == 0)
+ AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n");
+
}
//____________________________________________________________________________
// 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());
-
- AliPHOSGeometry * geom = gime->PHOSGeometry();
-
- fEmcCrystals = geom->GetNModules() * geom->GetNCristalsInModule() ;
+ fEmcCrystals = fGeom->GetNModules() * fGeom->GetNCristalsInModule() ;
if(!gMinuit)
gMinuit = new TMinuit(100);
- if ( !gime->Clusterizer() ) {
- gime->PostClusterizer(this);
- }
}
//____________________________________________________________________________
fNumberOfCpvClusters = 0 ;
fNumberOfEmcClusters = 0 ;
+
+ const AliPHOSRecoParam* parEmc = AliPHOSReconstructor::GetRecoParamEmc();
+ if(!parEmc) AliFatal("Reconstruction parameters for EMC not set!");
+
+ const AliPHOSRecoParam* parCpv = AliPHOSReconstructor::GetRecoParamCpv();
+ if(!parCpv) AliFatal("Reconstruction parameters for CPV not set!");
+
+ fCpvClusteringThreshold = parCpv->GetClusteringThreshold();
+ fEmcClusteringThreshold = parEmc->GetClusteringThreshold();
- fCpvClusteringThreshold = 0.0;
- fEmcClusteringThreshold = 0.2;
-
- fEmcLocMaxCut = 0.03 ;
- fCpvLocMaxCut = 0.03 ;
+ fEmcLocMaxCut = parEmc->GetLocalMaxCut();
+ fCpvLocMaxCut = parCpv->GetLocalMaxCut();
- fEmcMinE = 0.01 ;
- fCpvMinE = 0.0 ;
+ fEmcMinE = parEmc->GetMinE();
+ fCpvMinE = parCpv->GetMinE();
- fW0 = 4.5 ;
- fW0CPV = 4.0 ;
+ fW0 = parEmc->GetLogWeight();
+ fW0CPV = parCpv->GetLogWeight();
fEmcTimeGate = 1.e-8 ;
fToUnfold = kTRUE ;
- fRecPointsInRun = 0 ;
-
fWrite = kTRUE ;
fCalibData = 0 ;
- SetEventRange(0,-1) ;
-
fIsOldRCUFormat = kFALSE;
}
// 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] ;
- geom->AbsToRelNumbering(d1->GetId(), relid1) ;
+ fGeom->AbsToRelNumbering(d1->GetId(), relid1) ;
Int_t relid2[4] ;
- geom->AbsToRelNumbering(d2->GetId(), relid2) ;
+ fGeom->AbsToRelNumbering(d2->GetId(), relid2) ;
if ( (relid1[0] == relid2[0]) && (relid1[1]==relid2[1]) ) { // inside the same PHOS module
Int_t rowdiff = TMath::Abs( relid1[2] - relid2[2] ) ;
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 nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
+ Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
if(digit->GetId() <= nEMC ) rv = kTRUE;
Bool_t rv = kFALSE ;
- AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
-
- Int_t nEMC = geom->GetNModules()*geom->GetNPhi()*geom->GetNZ();
+ Int_t nEMC = fGeom->GetNModules()*fGeom->GetNPhi()*fGeom->GetNZ();
if(digit->GetId() > nEMC ) rv = kTRUE;
return rv ;
}
-//____________________________________________________________________________
-void AliPHOSClusterizerv1::Unload()
-{
- AliPHOSGetter * gime = AliPHOSGetter::Instance() ;
- gime->PhosLoader()->UnloadDigits() ;
- gime->PhosLoader()->UnloadRecPoints() ;
-}
-
//____________________________________________________________________________
void AliPHOSClusterizerv1::WriteRecPoints()
{
// Creates new branches with given title
// fills and writes into TreeR.
- AliPHOSGetter * gime = AliPHOSGetter::Instance();
-
- 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();
+ Int_t nEmc = fEMCRecPoints->GetEntriesFast();
for(index = 0; index < nEmc; index++){
- AliPHOSEmcRecPoint * rp = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) );
+ AliPHOSEmcRecPoint * rp =
+ dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) );
rp->Purify(fEmcMinE) ;
- if(rp->GetMultiplicity()>0.) //If this RP is not empty
- rp->EvalAll(fW0,digits) ;
- else{
- emcRecPoints->RemoveAt(index) ;
+ if(rp->GetMultiplicity()==0){
+ fEMCRecPoints->RemoveAt(index) ;
delete rp ;
+ continue;
}
+
+ // No vertex is available now, calculate corrections in PID
+ rp->EvalAll(fW0,fDigitsArr) ;
+ TVector3 fakeVtx(0.,0.,0.) ;
+ rp->EvalAll(fW0,fakeVtx,fDigitsArr) ;
+ rp->EvalLocal2TrackingCSTransform();
}
- emcRecPoints->Compress() ;
- emcRecPoints->Sort() ;
- // emcRecPoints->Expand(emcRecPoints->GetEntriesFast()) ;
- for(index = 0; index < emcRecPoints->GetEntries(); index++){
- dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) )->SetIndexInList(index) ;
+ fEMCRecPoints->Compress() ;
+ // fEMCRecPoints->Sort() ; //Can not sort until position is calculated!
+ // fEMCRecPoints->Expand(fEMCRecPoints->GetEntriesFast()) ;
+ for(index = 0; index < fEMCRecPoints->GetEntries(); index++){
+ dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) )->SetIndexInList(index) ;
}
+ //For each rec.point set the distance to the nearest bad crystal (BVP)
+ SetDistancesToBadChannels();
+
//Now the same for CPV
- for(index = 0; index < cpvRecPoints->GetEntries(); index++){
- AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) );
- rp->EvalAll(fW0CPV,digits) ;
+ for(index = 0; index < fCPVRecPoints->GetEntries(); index++){
+ AliPHOSCpvRecPoint * rp = dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) );
+ rp->EvalAll(fW0CPV,fDigitsArr) ;
+ rp->EvalLocal2TrackingCSTransform();
}
- cpvRecPoints->Sort() ;
+// fCPVRecPoints->Sort() ;
- for(index = 0; index < cpvRecPoints->GetEntries(); index++)
- dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(index) )->SetIndexInList(index) ;
+ for(index = 0; index < fCPVRecPoints->GetEntries(); index++)
+ dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(index) )->SetIndexInList(index) ;
- cpvRecPoints->Expand(cpvRecPoints->GetEntriesFast()) ;
+ fCPVRecPoints->Expand(fCPVRecPoints->GetEntriesFast()) ;
if(fWrite){ //We write TreeR
- TTree * treeR = gime->TreeR();
-
- Int_t bufferSize = 32000 ;
- Int_t splitlevel = 0 ;
-
- //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");
+ fTreeR->Fill();
}
}
// Steering method to construct the clusters stored in a list of Reconstructed Points
// A cluster is defined as a list of neighbour digits
-
- 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) ;
+ CleanDigits(fDigitsArr) ;
- TClonesArray * digitsC = static_cast<TClonesArray*>( digits->Clone() ) ;
+ TClonesArray * digitsC = static_cast<TClonesArray*>( fDigitsArr->Clone() ) ;
// Clusterization starts
if ( IsInEmc(digit) ) {
// start a new EMC RecPoint
- if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
- emcRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
+ if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
+ fEMCRecPoints->Expand(2*fNumberOfEmcClusters+1) ;
- emcRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
- clu = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ;
+ fEMCRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ;
+ clu = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
fNumberOfEmcClusters++ ;
clu->AddDigit(*digit, CalibrateEMC(digit->GetEnergy(),digit->GetId())) ;
clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
} else {
// start a new CPV cluster
- if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
- cpvRecPoints->Expand(2*fNumberOfCpvClusters+1);
+ if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
+ fCPVRecPoints->Expand(2*fNumberOfCpvClusters+1);
- cpvRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
+ fCPVRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ;
- clu = dynamic_cast<AliPHOSCpvRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;
+ clu = dynamic_cast<AliPHOSCpvRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
fNumberOfCpvClusters++ ;
clu->AddDigit(*digit, CalibrateCPV(digit->GetAmp(),digit->GetId()) ) ;
clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ;
AliPHOSDigit * digitN ;
index = 0 ;
while (index < iDigitInCluster){ // scan over digits already in cluster
- digit = dynamic_cast<AliPHOSDigit*>( digits->At(clusterdigitslist[index]) ) ;
+ digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->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 !!!!!
// Unfolds clusters using the shape of an ElectroMagnetic shower
// 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 = geom->GetNModules() ;
+ Int_t nModulesToUnfold = fGeom->GetNModules() ;
Int_t numberofNotUnfolded = fNumberOfEmcClusters ;
Int_t index ;
for(index = 0 ; index < numberofNotUnfolded ; index++){
- AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(index) ) ;
+ AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(index) ) ;
if(emcRecPoint->GetPHOSMod()> nModulesToUnfold)
break ;
Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
Float_t * maxAtEnergy = new Float_t[nMultipl] ;
- Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,digits) ;
+ Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fEmcLocMaxCut,fDigitsArr) ;
if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
- emcRecPoints->Remove(emcRecPoint);
- emcRecPoints->Compress() ;
+ fEMCRecPoints->Remove(emcRecPoint);
+ fEMCRecPoints->Compress() ;
index-- ;
fNumberOfEmcClusters -- ;
numberofNotUnfolded-- ;
// Unfold now CPV clusters
if(fNumberOfCpvClusters > 0){
- Int_t nModulesToUnfold = geom->GetNModules() ;
+ Int_t nModulesToUnfold = fGeom->GetNModules() ;
Int_t numberofCpvNotUnfolded = fNumberOfCpvClusters ;
Int_t index ;
for(index = 0 ; index < numberofCpvNotUnfolded ; index++){
- AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( cpvRecPoints->At(index) ) ;
+ AliPHOSRecPoint * recPoint = dynamic_cast<AliPHOSRecPoint *>( fCPVRecPoints->At(index) ) ;
if(recPoint->GetPHOSMod()> nModulesToUnfold)
break ;
Int_t nMultipl = emcRecPoint->GetMultiplicity() ;
AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ;
Float_t * maxAtEnergy = new Float_t[nMultipl] ;
- Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,digits) ;
+ Int_t nMax = emcRecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fCpvLocMaxCut,fDigitsArr) ;
if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0
UnfoldCluster(emcRecPoint, nMax, maxAt, maxAtEnergy) ;
- cpvRecPoints->Remove(emcRecPoint);
- cpvRecPoints->Compress() ;
+ fCPVRecPoints->Remove(emcRecPoint);
+ fCPVRecPoints->Compress() ;
index-- ;
numberofCpvNotUnfolded-- ;
fNumberOfCpvClusters-- ;
}
//____________________________________________________________________________
-Double_t AliPHOSClusterizerv1::ShowerShape(Double_t r)
+Double_t AliPHOSClusterizerv1::ShowerShape(Double_t x, Double_t z)
{
// Shape of the shower (see PHOS TDR)
// If you change this function, change also the gradient evaluation in ChiSquare()
- Double_t r4 = r*r*r*r ;
- Double_t r295 = TMath::Power(r, 2.95) ;
+ //for the moment we neglect dependence on the incident angle.
+
+ Double_t r2 = x*x + z*z ;
+ Double_t r4 = r2*r2 ;
+ Double_t r295 = TMath::Power(r2, 2.95/2.) ;
Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ;
return shape ;
}
{
// 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] ;
Int_t nDigits = iniEmc->GetMultiplicity() ;
Float_t * efit = new Float_t[nDigits] ;
- Float_t xDigit=0.,zDigit=0.,distance=0. ;
+ Float_t xDigit=0.,zDigit=0. ;
Float_t xpar=0.,zpar=0.,epar=0. ;
Int_t relid[4] ;
AliPHOSDigit * digit = 0 ;
Int_t * emcDigits = iniEmc->GetDigitsList() ;
+ TVector3 vIncid ;
+
Int_t iparam ;
Int_t iDigit ;
for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
- digit = dynamic_cast<AliPHOSDigit*>( digits->At(emcDigits[iDigit] ) ) ;
- geom->AbsToRelNumbering(digit->GetId(), relid) ;
- geom->RelPosInModule(relid, xDigit, zDigit) ;
+ digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At(emcDigits[iDigit] ) ) ;
+ fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ fGeom->RelPosInModule(relid, xDigit, zDigit) ;
efit[iDigit] = 0;
iparam = 0 ;
zpar = fitparameters[iparam+1] ;
epar = fitparameters[iparam+2] ;
iparam += 3 ;
- distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ;
- distance = TMath::Sqrt(distance) ;
- efit[iDigit] += epar * ShowerShape(distance) ;
+// fGeom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ;
+// efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) ;
+ efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
}
}
zpar = fitparameters[iparam+1] ;
epar = fitparameters[iparam+2] ;
iparam += 3 ;
+// fGeom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ;
AliPHOSEmcRecPoint * emcRP = 0 ;
- if(iniEmc->IsEmc()){ //create new entries in fEmcRecPoints...
+ if(iniEmc->IsEmc()){ //create new entries in fEMCRecPoints...
- if(fNumberOfEmcClusters >= emcRecPoints->GetSize())
- emcRecPoints->Expand(2*fNumberOfEmcClusters) ;
+ if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize())
+ fEMCRecPoints->Expand(2*fNumberOfEmcClusters) ;
- (*emcRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
- emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( emcRecPoints->At(fNumberOfEmcClusters) ) ;
+ (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ;
+ emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fEMCRecPoints->At(fNumberOfEmcClusters) ) ;
fNumberOfEmcClusters++ ;
emcRP->SetNExMax((Int_t)nPar/3) ;
}
- else{//create new entries in fCpvRecPoints
- if(fNumberOfCpvClusters >= cpvRecPoints->GetSize())
- cpvRecPoints->Expand(2*fNumberOfCpvClusters) ;
+ else{//create new entries in fCPVRecPoints
+ if(fNumberOfCpvClusters >= fCPVRecPoints->GetSize())
+ fCPVRecPoints->Expand(2*fNumberOfCpvClusters) ;
- (*cpvRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
- emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( cpvRecPoints->At(fNumberOfCpvClusters) ) ;
+ (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ;
+ emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( fCPVRecPoints->At(fNumberOfCpvClusters) ) ;
fNumberOfCpvClusters++ ;
}
Float_t eDigit ;
for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){
- 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] ;
+ digit = dynamic_cast<AliPHOSDigit*>( fDigitsArr->At( emcDigits[iDigit] ) ) ;
+ fGeom->AbsToRelNumbering(digit->GetId(), relid) ;
+ fGeom->RelPosInModule(relid, xDigit, zDigit) ;
+// ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ;
+ ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ;
eDigit = emcEnergies[iDigit] * ratio ;
emcRP->AddDigit( *digit, eDigit ) ;
}
AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint*>( toMinuit->At(0) ) ;
TClonesArray * digits = dynamic_cast<TClonesArray*>( toMinuit->At(1) ) ;
+ // A bit buggy way to get an access to the geometry
+ // To be revised!
+ AliPHOSGeometry *geom = dynamic_cast<AliPHOSGeometry *>(toMinuit->At(2));
-
+// TVector3 * vtx = dynamic_cast<TVector3*>(toMinuit->At(3)) ; //Vertex position
// AliPHOSEmcRecPoint * emcRP = dynamic_cast<AliPHOSEmcRecPoint *>( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit
Float_t * emcEnergies = emcRP->GetEnergiesList() ;
- const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ;
+// TVector3 vInc ;
fret = 0. ;
Int_t iparam ;
Int_t iParam = 0 ;
efit = 0 ;
while(iParam < nPar ){
- Double_t distance = (xDigit - x[iParam]) * (xDigit - x[iParam]) ;
+ Double_t dx = (xDigit - x[iParam]) ;
iParam++ ;
- distance += (zDigit - x[iParam]) * (zDigit - x[iParam]) ;
- distance = TMath::Sqrt( distance ) ;
+ Double_t dz = (zDigit - x[iParam]) ;
iParam++ ;
- efit += x[iParam] * ShowerShape(distance) ;
+// fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ;
+// efit += x[iParam] * ShowerShape(dx,dz,vInc) ;
+ efit += x[iParam] * ShowerShape(dx,dz) ;
iParam++ ;
}
Double_t sum = 2. * (efit - emcEnergies[iDigit]) / emcEnergies[iDigit] ; // Here we assume, that sigma = sqrt(E)
Double_t xpar = x[iParam] ;
Double_t zpar = x[iParam+1] ;
Double_t epar = x[iParam+2] ;
+// fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) );
- Double_t shape = sum * ShowerShape(dr) ;
+// Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
+ Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ;
+//DP: No incident angle dependence in gradient yet!!!!!!
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) ) +
Double_t zpar = x[iparam+1] ;
Double_t epar = x[iparam+2] ;
iparam += 3 ;
- Double_t distance = (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ;
- distance = TMath::Sqrt(distance) ;
- efit += epar * ShowerShape(distance) ;
+// fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ;
+// efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ;
+ efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ;
}
fret += (efit-emcEnergies[iDigit])*(efit-emcEnergies[iDigit])/emcEnergies[iDigit] ;
fCpvLocMaxCut,
fW0CPV )) ;
}
-
-
+//____________________________________________________________________________
+//void AliPHOSClusterizerv1::GetVertex(void)
+//{ //Extracts vertex posisition
+//
+ //ESD
+//DP - todo if(){
+//
+// }
+
+// //MC Generator
+// if(gAlice && gAlice->GetMCApp() && gAlice->Generator()){
+// Float_t x,y,z ;
+// gAlice->Generator()->GetOrigin(x,y,z) ;
+// fVtx.SetXYZ(x,y,z) ;
+// return ;
+// }
+//
+// //No any source
+// fVtx[0]=fVtx[1]=fVtx[2]=0. ;
+//
+//}
//____________________________________________________________________________
void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option)
{
// Prints list of RecPoints produced at the current pass of AliPHOSClusterizer
- 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() )) ;
+ AliInfo(Form("\nFound %d EMC RecPoints and %d CPV RecPoints",
+ fEMCRecPoints->GetEntriesFast(),
+ fCPVRecPoints->GetEntriesFast() )) ;
- fRecPointsInRun += emcRecPoints->GetEntriesFast() ;
- fRecPointsInRun += cpvRecPoints->GetEntriesFast() ;
-
-
if(strstr(option,"all")) {
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 < emcRecPoints->GetEntries() ; index++) {
- AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )emcRecPoints->At(index) ;
+ for (index = 0 ; index < fEMCRecPoints->GetEntries() ; index++) {
+ AliPHOSEmcRecPoint * rp = (AliPHOSEmcRecPoint * )fEMCRecPoints->At(index) ;
TVector3 locpos;
rp->GetLocalPosition(locpos);
Float_t lambda[2];
//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) ;
+ for (index = 0 ; index < fCPVRecPoints->GetEntries() ; index++) {
+ AliPHOSCpvRecPoint * rp = (AliPHOSCpvRecPoint * )fCPVRecPoints->At(index) ;
TVector3 locpos;
rp->GetLocalPosition(locpos);
}
}
+
+//____________________________________________________________________________
+void AliPHOSClusterizerv1::SetDistancesToBadChannels()
+{
+ //For each EMC rec. point set the distance to the nearest bad crystal.
+ //Author: Boris Polichtchouk
+
+ if(!fCalibData->GetNumOfEmcBadChannels()) return;
+ AliInfo(Form("%d bad channel(s) found.\n",fCalibData->GetNumOfEmcBadChannels()));
+
+ Int_t badIds[8000];
+ fCalibData->EmcBadChannelIds(badIds);
+
+ AliPHOSEmcRecPoint* rp;
+
+ TMatrixF gmat;
+ TVector3 gposRecPoint; // global (in ALICE frame) position of rec. point
+ TVector3 gposBadChannel; // global position of bad crystal
+ TVector3 dR;
+
+ Float_t dist,minDist;
+
+ for(Int_t iRP=0; iRP<fEMCRecPoints->GetEntries(); iRP++){
+ rp = (AliPHOSEmcRecPoint*)fEMCRecPoints->At(iRP);
+ minDist = 1.e+07;
+
+ for(Int_t iBad=0; iBad<fCalibData->GetNumOfEmcBadChannels(); iBad++) {
+ rp->GetGlobalPosition(gposRecPoint,gmat);
+ fGeom->RelPosInAlice(badIds[iBad],gposBadChannel);
+ AliDebug(2,Form("BC position:[%.3f,%.3f,%.3f], RP position:[%.3f,%.3f,%.3f]. E=%.3f\n",
+ gposBadChannel.X(),gposBadChannel.Y(),gposBadChannel.Z(),
+ gposRecPoint.X(),gposRecPoint.Y(),gposRecPoint.Z(),rp->GetEnergy()));
+ dR = gposBadChannel-gposRecPoint;
+ dist = dR.Mag();
+ if(dist<minDist) minDist = dist;
+ }
+
+ rp->SetDistanceToBadCrystal(minDist);
+ }
+
+}
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff