X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAliPHOSClusterizerv1.cxx;h=840ffd05454774e42831d4c4172e7c48c6b2a4ed;hb=d3b046c5c8af22f7e709ce6a32058e6726bb083c;hp=db6a36033bee465ba5e7c4a4b03cf4383e8085ef;hpb=cd79ec76c7b78c38e985468025ec48169d1923ce;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AliPHOSClusterizerv1.cxx b/PHOS/AliPHOSClusterizerv1.cxx index db6a36033be..840ffd05454 100644 --- a/PHOS/AliPHOSClusterizerv1.cxx +++ b/PHOS/AliPHOSClusterizerv1.cxx @@ -17,7 +17,41 @@ /* History of cvs commits: * - * $Log$ + * $Log: AliPHOSClusterizerv1.cxx,v $ + * Revision 1.118 2007/12/11 21:23:26 kharlov + * Added possibility to swith off unfolding + * + * Revision 1.117 2007/10/18 08:42:05 kharlov + * Bad channels cleaned before clusterization + * + * Revision 1.116 2007/10/01 20:24:08 kharlov + * Memory leaks fixed + * + * Revision 1.115 2007/09/26 14:22:17 cvetan + * Important changes to the reconstructor classes. Complete elimination of the run-loaders, which are now steered only from AliReconstruction. Removal of the corresponding Reconstruct() and FillESD() methods. + * + * Revision 1.114 2007/09/06 16:06:44 kharlov + * Absence of sorting results in loose of all unfolded clusters + * + * Revision 1.113 2007/08/28 12:55:07 policheh + * Loaders removed from the reconstruction code (C.Cheshkov) + * + * 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 * @@ -117,23 +151,15 @@ // 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 : 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() +// 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 --- @@ -141,24 +167,25 @@ #include "TMinuit.h" #include "TTree.h" #include "TBenchmark.h" +#include "TClonesArray.h" // --- Standard library --- // --- AliRoot header files --- #include "AliLog.h" -#include "AliRunLoader.h" -#include "AliGenerator.h" -#include "AliPHOSGetter.h" +#include "AliConfig.h" #include "AliPHOSGeometry.h" #include "AliPHOSClusterizerv1.h" #include "AliPHOSEmcRecPoint.h" #include "AliPHOSCpvRecPoint.h" #include "AliPHOSDigit.h" #include "AliPHOSDigitizer.h" -#include "AliPHOSCalibrationDB.h" #include "AliCDBManager.h" #include "AliCDBStorage.h" #include "AliCDBEntry.h" +#include "AliPHOSRecoParam.h" +#include "AliPHOSReconstructor.h" +#include "AliPHOSCalibData.h" ClassImp(AliPHOSClusterizerv1) @@ -166,53 +193,45 @@ ClassImp(AliPHOSClusterizerv1) AliPHOSClusterizerv1::AliPHOSClusterizerv1() : AliPHOSClusterizer(), 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), + fWrite(0), + fNumberOfEmcClusters(0), fNumberOfCpvClusters(0), + fEmcClusteringThreshold(0), fCpvClusteringThreshold(0), fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0), - fW0CPV(0), fRecPointsInRun(0), fEmcTimeGate(0), - fIsOldRCUFormat(0) + fW0CPV(0), + fTimeGateLowAmp(0.), fTimeGateLow(0.), fTimeGateHigh(0.), + fEcoreRadius(0) { // default ctor (to be used mainly by Streamer) - InitParameters() ; - fDefaultInit = kTRUE ; + fDefaultInit = kTRUE ; + + for(Int_t i=0; i<53760; i++){ + fDigitsUsed[i]=0 ; + } } //____________________________________________________________________________ -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), + fWrite(0), + fNumberOfEmcClusters(0), fNumberOfCpvClusters(0), + fEmcClusteringThreshold(0), fCpvClusteringThreshold(0), fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0), - fW0CPV(0), fRecPointsInRun(0), fEmcTimeGate(0), - fIsOldRCUFormat(0) + fW0CPV(0), + fTimeGateLowAmp(0.), fTimeGateLow(0.), fTimeGateHigh(0.), + fEcoreRadius(0) { // ctor with the indication of the file where header Tree and digits Tree are stored - InitParameters() ; + for(Int_t i=0; i<53760; i++){ + fDigitsUsed[i]=0 ; + } + Init() ; 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() { @@ -220,71 +239,11 @@ AliPHOSClusterizerv1::AliPHOSClusterizerv1(const AliPHOSClusterizerv1 & obj) : } //____________________________________________________________________________ -const TString AliPHOSClusterizerv1::BranchName() const -{ - return GetName(); -} - -//____________________________________________________________________________ -Float_t AliPHOSClusterizerv1::CalibrateEMC(Float_t amp, Int_t absId) -{ - // 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; -} - -//____________________________________________________________________________ -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"); @@ -294,60 +253,28 @@ void AliPHOSClusterizerv1::Exec(Option_t *option) return ; } - 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 { - AliRunLoader * rl = AliRunLoader::GetRunLoader(gime->PhosLoader()->GetTitle()); - rl->GetEvent(ievent); - 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() ; - - WriteRecPoints(); + if(fToUnfold) + MakeUnfolding(); - if(strstr(option,"deb")) - PrintRecPoints(option) ; + 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"); - 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"))); + } + fEMCRecPoints->Delete(); + fCPVRecPoints->Delete(); } //____________________________________________________________________________ @@ -359,9 +286,8 @@ Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers - AliPHOSGetter * gime = AliPHOSGetter::Instance(); - TClonesArray * digits = gime->Digits(); - + if(!gMinuit) //it was deleted by someone else + gMinuit = new TMinuit(100) ; gMinuit->mncler(); // Reset Minuit's list of paramters gMinuit->SetPrintLevel(-1) ; // No Printout gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ; @@ -369,7 +295,8 @@ Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** 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 @@ -382,16 +309,14 @@ Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** 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] ; @@ -443,23 +368,6 @@ Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** } -//____________________________________________________________________________ -void AliPHOSClusterizerv1::GetCalibrationParameters() -{ - // Set calibration parameters: - // if calibration database exists, they are read from database, - // 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"); - - fCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number - 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"); - -} //____________________________________________________________________________ void AliPHOSClusterizerv1::Init() @@ -467,20 +375,18 @@ 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()); - - 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); - } + if (!fgCalibData) + fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number + if (fgCalibData->GetCalibDataEmc() == 0) + AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n"); + if (fgCalibData->GetCalibDataCpv() == 0) + AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n"); + } //____________________________________________________________________________ @@ -489,32 +395,30 @@ 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 ; + const AliPHOSRecoParam* recoParam = AliPHOSReconstructor::GetRecoParam(); + if(!recoParam) AliFatal("Reconstruction parameters are not set!"); - fW0 = 4.5 ; - fW0CPV = 4.0 ; + recoParam->Print(); - fEmcTimeGate = 1.e-8 ; + fEmcClusteringThreshold = recoParam->GetEMCClusteringThreshold(); + fCpvClusteringThreshold = recoParam->GetCPVClusteringThreshold(); - fToUnfold = kTRUE ; - - fRecPointsInRun = 0 ; - - fWrite = kTRUE ; + fEmcLocMaxCut = recoParam->GetEMCLocalMaxCut(); + fCpvLocMaxCut = recoParam->GetCPVLocalMaxCut(); - fCalibData = 0 ; + fW0 = recoParam->GetEMCLogWeight(); + fW0CPV = recoParam->GetCPVLogWeight(); - SetEventRange(0,-1) ; + fTimeGateLowAmp = recoParam->GetTimeGateAmpThresh() ; + fTimeGateLow = recoParam->GetTimeGateLow() ; + fTimeGateHigh = recoParam->GetTimeGateHigh() ; - fIsOldRCUFormat = kFALSE; + fEcoreRadius = recoParam->GetEMCEcoreRadius(); + + fToUnfold = recoParam->EMCToUnfold() ; + + fWrite = kTRUE ; } //____________________________________________________________________________ @@ -523,67 +427,56 @@ Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)c // Gives the neighbourness of two digits = 0 are not neighbour but continue searching // = 1 are neighbour // = 2 are not neighbour but do not continue searching + // =-1 are not neighbour, continue searching, but do not look before d2 next time // neighbours are defined as digits having at least a common vertex // 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] ) ; 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 ; + if (( coldiff <= 1 ) && ( rowdiff <= 1 )){ //At least common vertex + // if (( relid1[2]==relid2[2] && coldiff <= 1 ) || ( relid1[3]==relid2[3] && rowdiff <= 1 )){ //common side + if((relid1[1] != 0) || CheckTimeGate(d1->GetTime(),d1->GetEnergy(),d2->GetTime(),d2->GetEnergy())) + return 1 ; } else { if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1)) - rv = 2; // Difference in row numbers is too large to look further + return 2; // Difference in row numbers is too large to look further } + return 0 ; } else { + if(relid1[0] > relid2[0] && relid1[1]==relid2[1] ) //we switched to the next module + return -1 ; + if(relid1[1] < relid2[1]) //we switched from EMC(0) to CPV(-1) + return -1 ; - if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) ) - rv=2 ; + return 2 ; } - return rv ; + return 0 ; } //____________________________________________________________________________ -void AliPHOSClusterizerv1::CleanDigits(TClonesArray * digits) -{ - // Remove digits with amplitudes below threshold - - for(Int_t i=0; iGetEntriesFast(); i++){ - AliPHOSDigit * digit = static_cast(digits->At(i)) ; - if ( (IsInEmc(digit) && CalibrateEMC(digit->GetEnergy(),digit->GetId()) < fEmcMinE) || - (IsInCpv(digit) && CalibrateCPV(digit->GetAmp() ,digit->GetId()) < fCpvMinE) ) - digits->RemoveAt(i) ; +Bool_t AliPHOSClusterizerv1::CheckTimeGate(Float_t t1, Float_t amp1, Float_t t2, Float_t amp2)const{ + //Check if two cells have reasonable time difference + //Note that at low amplitude time is defined up to 1 tick == 100 ns. + if(amp1Compress() ; - for (Int_t i = 0 ; i < digits->GetEntriesFast() ; i++) { - AliPHOSDigit *digit = static_cast( digits->At(i) ) ; - digit->SetIndexInList(i) ; + else{ //Time should be measured with good accuracy + return (TMath::Abs(t1 - t2 ) < fTimeGateHigh) ; } - //Overwrite digits tree - AliPHOSGetter* gime = AliPHOSGetter::Instance(); - TTree * treeD = gime->TreeD(); - treeD->Branch("PHOS", &digits); - treeD->Fill() ; - gime->WriteDigits("OVERWRITE"); - gime->PhosLoader()->UnloadDigits() ; } //____________________________________________________________________________ Bool_t AliPHOSClusterizerv1::IsInEmc(AliPHOSDigit * digit) const @@ -591,9 +484,8 @@ 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; @@ -607,23 +499,13 @@ Bool_t AliPHOSClusterizerv1::IsInCpv(AliPHOSDigit * digit) const 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() { @@ -631,73 +513,53 @@ 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(); + Float_t emcMinE= AliPHOSReconstructor::GetRecoParam()->GetEMCMinE(); //Minimal digit energy + TVector3 fakeVtx(0.,0.,0.) ; for(index = 0; index < nEmc; index++){ AliPHOSEmcRecPoint * rp = - dynamic_cast( emcRecPoints->At(index) ); - rp->Purify(fEmcMinE) ; + static_cast( fEMCRecPoints->At(index) ); + rp->Purify(emcMinE,fDigitsArr) ; if(rp->GetMultiplicity()==0){ - emcRecPoints->RemoveAt(index) ; + fEMCRecPoints->RemoveAt(index) ; delete rp ; continue; } // No vertex is available now, calculate corrections in PID - rp->EvalAll(fW0,digits) ; - TVector3 fakeVtx(0.,0.,0.) ; - rp->EvalAll(fW0,fakeVtx,digits) ; + rp->EvalAll(fDigitsArr) ; + rp->EvalCoreEnergy(fW0,fEcoreRadius,fDigitsArr) ; + rp->EvalAll(fW0,fakeVtx,fDigitsArr) ; rp->EvalLocal2TrackingCSTransform(); } - emcRecPoints->Compress() ; -// emcRecPoints->Sort() ; //Can not sort until position is calculated! - // emcRecPoints->Expand(emcRecPoints->GetEntriesFast()) ; - for(index = 0; index < emcRecPoints->GetEntries(); index++){ - dynamic_cast( emcRecPoints->At(index) )->SetIndexInList(index) ; + fEMCRecPoints->Compress() ; + fEMCRecPoints->Sort() ; + // fEMCRecPoints->Expand(fEMCRecPoints->GetEntriesFast()) ; + for(index = 0; index < fEMCRecPoints->GetEntries(); index++){ + static_cast( 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( cpvRecPoints->At(index) ); - rp->EvalAll(fW0CPV,digits) ; + for(index = 0; index < fCPVRecPoints->GetEntries(); index++){ + AliPHOSCpvRecPoint * rp = static_cast( fCPVRecPoints->At(index) ); + rp->EvalAll(fDigitsArr) ; + rp->EvalAll(fW0CPV,fakeVtx,fDigitsArr) ; rp->EvalLocal2TrackingCSTransform(); } -// cpvRecPoints->Sort() ; + fCPVRecPoints->Sort() ; - for(index = 0; index < cpvRecPoints->GetEntries(); index++) - dynamic_cast( cpvRecPoints->At(index) )->SetIndexInList(index) ; + for(index = 0; index < fCPVRecPoints->GetEntries(); index++) + static_cast( 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(); } } @@ -707,124 +569,101 @@ void AliPHOSClusterizerv1::MakeClusters() // 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) ; + fNumberOfCpvClusters = 0 ; + fNumberOfEmcClusters = 0 ; - TClonesArray * digitsC = static_cast( digits->Clone() ) ; - - // Clusterization starts + //Mark all digits as unused yet + const Int_t maxNDigits = 3584; // There is no clusters larger than PHOS module ;) + Int_t nDigits=fDigitsArr->GetEntriesFast() ; - TIter nextdigit(digitsC) ; + for(Int_t i=0; i( nextdigit()) ) ) { // scan over the list of digitsC - + digit=static_cast(fDigitsArr->At(i)) ; - AliPHOSRecPoint * clu = 0 ; + clu=0 ; - TArrayI clusterdigitslist(1500) ; Int_t index ; - if (( IsInEmc (digit) && - CalibrateEMC(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) || - ( IsInCpv (digit) && - CalibrateCPV(digit->GetAmp() ,digit->GetId()) > fCpvClusteringThreshold ) ) { + //is this digit so energetic that start cluster? + if (( IsInEmc(digit) && Calibrate(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) || + ( IsInCpv(digit) && Calibrate(digit->GetEnergy(),digit->GetId()) > fCpvClusteringThreshold ) ) { Int_t iDigitInCluster = 0 ; - 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( emcRecPoints->At(fNumberOfEmcClusters) ) ; + fEMCRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ; + clu = static_cast( fEMCRecPoints->At(fNumberOfEmcClusters) ) ; fNumberOfEmcClusters++ ; - clu->AddDigit(*digit, CalibrateEMC(digit->GetEnergy(),digit->GetId())) ; + clu->AddDigit(*digit, Calibrate(digit->GetEnergy(),digit->GetId()),CalibrateT(digit->GetTime(),digit->GetId(),digit->IsLG())) ; clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ; iDigitInCluster++ ; - digitsC->Remove(digit) ; - + fDigitsUsed[i]=kTRUE ; } 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) ; - - clu = dynamic_cast( cpvRecPoints->At(fNumberOfCpvClusters) ) ; + fCPVRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ; + clu = static_cast( fCPVRecPoints->At(fNumberOfCpvClusters) ) ; fNumberOfCpvClusters++ ; - clu->AddDigit(*digit, CalibrateCPV(digit->GetAmp(),digit->GetId()) ) ; + clu->AddDigit(*digit, Calibrate(digit->GetEnergy(),digit->GetId()),0.) ; // no timing information in CPV clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ; iDigitInCluster++ ; - digitsC->Remove(digit) ; - nextdigit.Reset() ; - - // Here we remove remaining EMC digits, which cannot make a cluster - - if( notremoved ) { - while( ( digit = dynamic_cast( nextdigit() ) ) ) { - if( IsInEmc(digit) ) - digitsC->Remove(digit) ; - else - break ; - } - notremoved = kFALSE ; - } - + fDigitsUsed[i]=kTRUE ; } // else - - nextdigit.Reset() ; - + + //Now scan remaining digits in list to find neigbours of our seed + AliPHOSDigit * digitN ; index = 0 ; while (index < iDigitInCluster){ // scan over digits already in cluster - digit = dynamic_cast( digits->At(clusterdigitslist[index]) ) ; + digit = static_cast( fDigitsArr->At(clusterdigitslist[index]) ) ; index++ ; - while ( (digitN = dynamic_cast( nextdigit() ) ) ) { // scan over the reduced list of digits + for(Int_t j=iFirst; j= maxNDigits) { + AliError(Form("The number of digits in cluster is more than %d, skip the rest of event", + maxNDigits)); + return; + } + if(fDigitsUsed[j]) + continue ; //look through remaining digits + digitN = static_cast( fDigitsArr->At(j) ) ; Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!! switch (ineb ) { + case -1: //too early (e.g. previous module), do not look before j at subsequent passes + iFirst=j ; + break ; case 0 : // not a neighbour break ; case 1 : // are neighbours - if (IsInEmc (digitN)) - clu->AddDigit(*digitN, CalibrateEMC( digitN->GetEnergy(), digitN->GetId() ) ); - else - clu->AddDigit(*digitN, CalibrateCPV( digitN->GetAmp() , digitN->GetId() ) ); - - clusterdigitslist[iDigitInCluster] = digitN->GetIndexInList() ; + clu->AddDigit(*digitN, Calibrate(digitN->GetEnergy(),digitN->GetId()),CalibrateT(digitN->GetTime(),digitN->GetId(),digit->IsLG())) ; + clusterdigitslist[iDigitInCluster] = j ; iDigitInCluster++ ; - digitsC->Remove(digitN) ; + fDigitsUsed[j]=kTRUE ; break ; case 2 : // too far from each other - goto endofloop; + goto endOfLoop; } // switch - } // while digitN + } - endofloop: ; - nextdigit.Reset() ; + endOfLoop: ; //scanned all possible neighbours for this digit } // loop over cluster - } // energy theshold - - - } // while digit - - delete digitsC ; + } } @@ -834,36 +673,29 @@ void AliPHOSClusterizerv1::MakeUnfolding() // 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( emcRecPoints->At(index) ) ; + AliPHOSEmcRecPoint * emcRecPoint = static_cast( 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-- ; @@ -882,28 +714,28 @@ void AliPHOSClusterizerv1::MakeUnfolding() // 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( cpvRecPoints->At(index) ) ; + AliPHOSRecPoint * recPoint = static_cast( fCPVRecPoints->At(index) ) ; if(recPoint->GetPHOSMod()> nModulesToUnfold) break ; - AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast(recPoint) ; + AliPHOSEmcRecPoint * emcRecPoint = static_cast(recPoint) ; 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-- ; @@ -940,18 +772,11 @@ void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc, { // 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) ; @@ -976,9 +801,9 @@ void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc, Int_t iparam ; Int_t iDigit ; for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ - digit = dynamic_cast( digits->At(emcDigits[iDigit] ) ) ; - geom->AbsToRelNumbering(digit->GetId(), relid) ; - geom->RelPosInModule(relid, xDigit, zDigit) ; + digit = static_cast( fDigitsArr->At(emcDigits[iDigit] ) ) ; + fGeom->AbsToRelNumbering(digit->GetId(), relid) ; + fGeom->RelPosInModule(relid, xDigit, zDigit) ; efit[iDigit] = 0; iparam = 0 ; @@ -987,13 +812,12 @@ void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc, zpar = fitparameters[iparam+1] ; epar = fitparameters[iparam+2] ; iparam += 3 ; -// geom->GetIncidentVector(fVtx,relid[0],xpar,zpar,vIncid) ; +// 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) ; } } - // Now create new RecPoints and fill energy lists with efit corrected to fluctuations // so that energy deposited in each cell is distributed betwin new clusters proportionally // to its contribution to efit @@ -1007,38 +831,38 @@ void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc, zpar = fitparameters[iparam+1] ; epar = fitparameters[iparam+2] ; iparam += 3 ; -// geom->GetIncidentVector(fVtx,iniEmc->GetPHOSMod(),xpar,zpar,vIncid) ; - +// 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( emcRecPoints->At(fNumberOfEmcClusters) ) ; + (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ; + emcRP = static_cast( 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( cpvRecPoints->At(fNumberOfCpvClusters) ) ; + (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ; + emcRP = static_cast( fCPVRecPoints->At(fNumberOfCpvClusters) ) ; fNumberOfCpvClusters++ ; } Float_t eDigit ; for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ - digit = dynamic_cast( digits->At( emcDigits[iDigit] ) ) ; - geom->AbsToRelNumbering(digit->GetId(), relid) ; - geom->RelPosInModule(relid, xDigit, zDigit) ; + digit = static_cast( 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 ) ; + emcRP->AddDigit( *digit, eDigit,CalibrateT(digit->GetTime(),digit->GetId(),digit->IsLG()) ) ; } } @@ -1053,13 +877,17 @@ void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Dou // Calculates the Chi square for the cluster unfolding minimization // Number of parameters, Gradient, Chi squared, parameters, what to do - TList * toMinuit = dynamic_cast( gMinuit->GetObjectFit() ) ; + TList * toMinuit = static_cast( gMinuit->GetObjectFit() ) ; + + AliPHOSEmcRecPoint * emcRP = static_cast( toMinuit->At(0) ) ; + TClonesArray * digits = static_cast( toMinuit->At(1) ) ; + // A bit buggy way to get an access to the geometry + // To be revised! + AliPHOSGeometry *geom = static_cast(toMinuit->At(2)); - AliPHOSEmcRecPoint * emcRP = dynamic_cast( toMinuit->At(0) ) ; - TClonesArray * digits = dynamic_cast( toMinuit->At(1) ) ; -// TVector3 * vtx = dynamic_cast(toMinuit->At(2)) ; //Vertex position +// TVector3 * vtx = static_cast(toMinuit->At(3)) ; //Vertex position - // AliPHOSEmcRecPoint * emcRP = dynamic_cast( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit + // AliPHOSEmcRecPoint * emcRP = static_cast( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit Int_t * emcDigits = emcRP->GetDigitsList() ; @@ -1067,7 +895,6 @@ void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Dou Float_t * emcEnergies = emcRP->GetEnergiesList() ; - const AliPHOSGeometry * geom = AliPHOSGetter::Instance()->PHOSGeometry() ; // TVector3 vInc ; fret = 0. ; Int_t iparam ; @@ -1083,7 +910,7 @@ void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Dou for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) { - digit = dynamic_cast( digits->At( emcDigits[iDigit] ) ); + digit = static_cast( digits->At( emcDigits[iDigit] ) ); Int_t relid[4] ; Float_t xDigit ; @@ -1101,7 +928,7 @@ void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Dou iParam++ ; Double_t dz = (zDigit - x[iParam]) ; iParam++ ; -// geom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),x[iParam-2],x[iParam-1],vInc) ; +// 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++ ; @@ -1112,7 +939,7 @@ void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Dou Double_t xpar = x[iParam] ; Double_t zpar = x[iParam+1] ; Double_t epar = x[iParam+2] ; -// geom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ; +// 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(xDigit - xpar,zDigit - zpar,vInc) ; Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ; @@ -1138,7 +965,7 @@ void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Dou Double_t zpar = x[iparam+1] ; Double_t epar = x[iparam+2] ; iparam += 3 ; -// geom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ; +// fGeom->GetIncidentVector(*vtx,emcRP->GetPHOSMod(),xpar,zpar,vInc) ; // efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar,vInc) ; efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ; } @@ -1179,7 +1006,7 @@ void AliPHOSClusterizerv1::Print(const Option_t *)const else message = " AliPHOSClusterizerv1 not initialized " ; - AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(), + AliInfo(Form("%s, %s %s %s %s %f %f %f %f %f %f", message.Data(), taskName.Data(), GetTitle(), taskName.Data(), @@ -1192,51 +1019,20 @@ void AliPHOSClusterizerv1::Print(const Option_t *)const 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]; @@ -1257,8 +1053,8 @@ void AliPHOSClusterizerv1::PrintRecPoints(Option_t * option) //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); @@ -1277,16 +1073,11 @@ 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())); + if(!fgCalibData->GetNumOfEmcBadChannels()) return; - AliPHOSGetter* gime = AliPHOSGetter::Instance(); - AliPHOSGeometry* geom = gime->PHOSGeometry(); - - TObjArray * emcRecPoints = gime->EmcRecPoints() ; - Int_t badIds[8000]; - fCalibData->EmcBadChannelIds(badIds); + memset(badIds,0,8000*sizeof(Int_t)); + fgCalibData->EmcBadChannelIds(badIds); AliPHOSEmcRecPoint* rp; @@ -1296,14 +1087,24 @@ void AliPHOSClusterizerv1::SetDistancesToBadChannels() TVector3 dR; Float_t dist,minDist; - - for(Int_t iRP=0; iRPGetEntries(); iRP++){ - rp = (AliPHOSEmcRecPoint*)emcRecPoints->At(iRP); - minDist = 1.e+07; - - for(Int_t iBad=0; iBadGetNumOfEmcBadChannels(); iBad++) { + Int_t relid[4]={0,0,0,0} ; + TVector3 lpos ; + for(Int_t iRP=0; iRPGetEntries(); iRP++){ + rp = (AliPHOSEmcRecPoint*)fEMCRecPoints->At(iRP); + //evaluate distance to border + relid[0]=rp->GetPHOSMod() ; + relid[2]=1 ; + relid[3]=1 ; + Float_t xcorner,zcorner; + fGeom->RelPosInModule(relid,xcorner,zcorner) ; //coordinate of the corner cell + rp->GetLocalPosition(lpos) ; + minDist = 2.2+TMath::Min(-xcorner-TMath::Abs(lpos.X()),-zcorner-TMath::Abs(lpos.Z())); //2.2 - crystal size + for(Int_t iBad=0; iBadGetNumOfEmcBadChannels(); iBad++) { + fGeom->AbsToRelNumbering(badIds[iBad],relid) ; + if(relid[0]!=rp->GetPHOSMod()) //We can not evaluate global position directly since + continue ; //bad channels can be in the module which does not exist in simulations. rp->GetGlobalPosition(gposRecPoint,gmat); - geom->RelPosInAlice(badIds[iBad],gposBadChannel); + 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())); @@ -1316,3 +1117,48 @@ void AliPHOSClusterizerv1::SetDistancesToBadChannels() } } +//================================================================================== +Float_t AliPHOSClusterizerv1::Calibrate(Float_t amp, Int_t absId) const{ + // Calibrate EMC digit, i.e. multiply its Amp by a factor read from CDB + + const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ; + + //Determine rel.position of the cell absolute ID + Int_t relId[4]; + geom->AbsToRelNumbering(absId,relId); + Int_t module=relId[0]; + Int_t row =relId[2]; + Int_t column=relId[3]; + if(relId[1]){ //CPV + Float_t calibration = fgCalibData->GetADCchannelCpv(module,column,row); + return amp*calibration ; + } + else{ //EMC + Float_t calibration = fgCalibData->GetADCchannelEmc(module,column,row); + return amp*calibration ; + } +} +//================================================================================== +Float_t AliPHOSClusterizerv1::CalibrateT(Float_t time, Int_t absId,Bool_t isLG)const{ + // Calibrate time in EMC digit + + const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ; + + //Determine rel.position of the cell absolute ID + Int_t relId[4]; + geom->AbsToRelNumbering(absId,relId); + Int_t module=relId[0]; + Int_t row =relId[2]; + Int_t column=relId[3]; + if(relId[1]){ //CPV + return 0. ; + } + else{ //EMC + if(isLG) + time += fgCalibData->GetLGTimeShiftEmc(module,column,row); + else + time += fgCalibData->GetTimeShiftEmc(module,column,row); + return time ; + } +} +