X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliEMCALClusterizerv1.cxx;h=9cd9bcdda9309a3a8a8b3c626e561ed76faf9660;hb=dd70110949cdcdc8c7c1cba590a9a2ea7e6de638;hp=dce0d2fa9c3b7bfbb2e6b5fa7d7c27ef5c1dfc69;hpb=a5c60732b9c666d5300eae697ebe76a26883d3ab;p=u%2Fmrichter%2FAliRoot.git diff --git a/EMCAL/AliEMCALClusterizerv1.cxx b/EMCAL/AliEMCALClusterizerv1.cxx index dce0d2fa9c3..9cd9bcdda93 100644 --- a/EMCAL/AliEMCALClusterizerv1.cxx +++ b/EMCAL/AliEMCALClusterizerv1.cxx @@ -18,6 +18,7 @@ //*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute) // August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction // of new IO (à la PHOS) +// Mar 2007, Aleksei Pavlinov - new algoritmh of pseudo clusters ////////////////////////////////////////////////////////////////////////////// // Clusterization class. Performs clusterization (collects neighbouring active cells) and // unfolds the clusters having several local maxima. @@ -45,23 +46,27 @@ // // time - print benchmarking results // --- ROOT system --- +#include -#include +class TROOT; #include #include -#include +class TFolder; #include #include #include -#include +class TSystem; #include #include +#include + // --- Standard library --- // --- AliRoot header files --- #include "AliRunLoader.h" #include "AliRun.h" +#include "AliESD.h" #include "AliEMCALLoader.h" #include "AliEMCALClusterizerv1.h" #include "AliEMCALRecPoint.h" @@ -69,33 +74,82 @@ #include "AliEMCALDigitizer.h" #include "AliEMCAL.h" #include "AliEMCALGeometry.h" +#include "AliEMCALRawUtils.h" #include "AliEMCALHistoUtilities.h" #include "AliCDBManager.h" -#include "AliCDBStorage.h" + +class AliCDBStorage; #include "AliCDBEntry.h" ClassImp(AliEMCALClusterizerv1) //____________________________________________________________________________ - AliEMCALClusterizerv1::AliEMCALClusterizerv1() : AliEMCALClusterizer() +AliEMCALClusterizerv1::AliEMCALClusterizerv1() + : AliEMCALClusterizer(), + fHists(0),fPointE(0),fPointL1(0),fPointL2(0), + fPointDis(0),fPointMult(0),fDigitAmp(0),fMaxE(0), + fMaxL1(0),fMaxL2(0),fMaxDis(0),fGeom(0), + fDefaultInit(kTRUE), + fToUnfold(kFALSE), + fNumberOfECAClusters(0),fNTowerInGroup(0),fCalibData(0), + fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.), + fECAW0(0.),fRecPointsInRun(0),fTimeCut(0.),fMinECut(0.) { // default ctor (to be used mainly by Streamer) InitParameters() ; - fDefaultInit = kTRUE ; - fGeom = AliEMCALGeometry::GetInstance(); + fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaulGeometryName()); fGeom->GetTransformationForSM(); // Global <-> Local } //____________________________________________________________________________ AliEMCALClusterizerv1::AliEMCALClusterizerv1(const TString alirunFileName, const TString eventFolderName) -:AliEMCALClusterizer(alirunFileName, eventFolderName) + : AliEMCALClusterizer(alirunFileName, eventFolderName), + fHists(0),fPointE(0),fPointL1(0),fPointL2(0), + fPointDis(0),fPointMult(0),fDigitAmp(0),fMaxE(0), + fMaxL1(0),fMaxL2(0),fMaxDis(0),fGeom(0), + fDefaultInit(kFALSE), + fToUnfold(kFALSE), + fNumberOfECAClusters(0),fNTowerInGroup(0),fCalibData(0), + fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.), + fECAW0(0.),fRecPointsInRun(0),fTimeCut(0.),fMinECut(0.) { // ctor with the indication of the file where header Tree and digits Tree are stored InitParameters() ; Init() ; - fDefaultInit = kFALSE ; +} + +//____________________________________________________________________________ +AliEMCALClusterizerv1::AliEMCALClusterizerv1(const AliEMCALClusterizerv1& clus) + : AliEMCALClusterizer(clus), + fHists(clus.fHists), + fPointE(clus.fPointE), + fPointL1(clus.fPointL1), + fPointL2(clus.fPointL2), + fPointDis(clus.fPointDis), + fPointMult(clus.fPointMult), + fDigitAmp(clus.fDigitAmp), + fMaxE(clus.fMaxE), + fMaxL1(clus.fMaxL1), + fMaxL2(clus.fMaxL2), + fMaxDis(clus.fMaxDis), + fGeom(clus.fGeom), + fDefaultInit(clus.fDefaultInit), + fToUnfold(clus.fToUnfold), + fNumberOfECAClusters(clus.fNumberOfECAClusters), + fNTowerInGroup(clus.fNTowerInGroup), + fCalibData(clus.fCalibData), + fADCchannelECA(clus.fADCchannelECA), + fADCpedestalECA(clus.fADCpedestalECA), + fECAClusteringThreshold(clus.fECAClusteringThreshold), + fECALocMaxCut(clus.fECALocMaxCut), + fECAW0(clus.fECAW0), + fRecPointsInRun(clus.fRecPointsInRun), + fTimeCut(clus.fTimeCut), + fMinECut(clus.fMinECut) +{ + //copy ctor } //____________________________________________________________________________ @@ -119,37 +173,30 @@ Float_t AliEMCALClusterizerv1::Calibrate(Int_t amp, Int_t AbsId) // or from digitizer parameters for simulated data. if(fCalibData){ - - //JLK 13-Mar-2006 - //We now get geometry at a higher level - // - // Loader - // AliRunLoader *rl = AliRunLoader::GetRunLoader(); - - // Load EMCAL Geomtry - // rl->LoadgAlice(); - //AliRun * gAlice = rl->GetAliRun(); - //AliEMCAL * emcal = (AliEMCAL*)gAlice->GetDetector("EMCAL"); - //AliEMCALGeometry * geom = emcal->GetGeometry(); if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader") ; Int_t iSupMod = -1; - Int_t nTower = -1; + Int_t nModule = -1; Int_t nIphi = -1; Int_t nIeta = -1; Int_t iphi = -1; Int_t ieta = -1; - Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nTower, nIphi, nIeta) ; - if(!bCell) - Error("DigitizeEnergy","Wrong cell id number") ; - fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nTower,nIphi, nIeta,iphi,ieta); + Bool_t bCell = fGeom->GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ; + if(!bCell) { + fGeom->PrintGeometry(); + Error("Calibrate()"," Wrong cell id number : %i", AbsId); + assert(0); + } + + fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta); fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi); fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi); - return -fADCpedestalECA + amp * fADCchannelECA ; + + return -fADCpedestalECA + amp * fADCchannelECA ; } else //Return energy with default parameters if calibration is not available @@ -160,59 +207,47 @@ Float_t AliEMCALClusterizerv1::Calibrate(Int_t amp, Int_t AbsId) //____________________________________________________________________________ void AliEMCALClusterizerv1::Exec(Option_t * option) { - // Steering method to perform clusterization for events - // in the range from fFirstEvent to fLastEvent. - // This range is optionally set by SetEventRange(). - // if fLastEvent=-1 (by default), then process events until the end. + // Steering method to perform clusterization for the current event + // in AliEMCALLoader if(strstr(option,"tim")) gBenchmark->Start("EMCALClusterizer"); if(strstr(option,"print")) Print("") ; - + AliRunLoader *rl = AliRunLoader::GetRunLoader(); AliEMCALLoader *emcalLoader = dynamic_cast(rl->GetDetectorLoader("EMCAL")); //Get calibration parameters from file or digitizer default values. GetCalibrationParameters() ; - if (fLastEvent == -1) - fLastEvent = rl->GetNumberOfEvents() - 1; - Int_t nEvents = fLastEvent - fFirstEvent + 1; - Int_t ievent ; - rl->LoadDigits("EMCAL"); - for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) { - rl->GetEvent(ievent); + fNumberOfECAClusters = 0; - fNumberOfECAClusters = 0; + if(strstr(option,"pseudo")) + MakeClusters("pseudo") ; //both types + else + MakeClusters("") ; //only the real clusters - if(strstr(option,"pseudo")) - MakeClusters("pseudo") ; //both types - else - MakeClusters("") ; //only the real clusters + if(fToUnfold) + MakeUnfolding() ; - if(fToUnfold) - MakeUnfolding() ; + WriteRecPoints() ; - WriteRecPoints() ; + if(strstr(option,"deb") || strstr(option,"all")) + PrintRecPoints(option) ; - if(strstr(option,"deb")) - PrintRecPoints(option) ; + AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",emcalLoader->RecPoints()->GetEntriesFast())); - //increment the total number of recpoints per run - fRecPointsInRun += emcalLoader->RecPoints()->GetEntriesFast() ; - } - - Unload(); + //increment the total number of recpoints per run + fRecPointsInRun += emcalLoader->RecPoints()->GetEntriesFast() ; if(strstr(option,"tim")){ gBenchmark->Stop("EMCALClusterizer"); - printf("Exec took %f seconds for Clusterizing %f seconds per event", - gBenchmark->GetCpuTime("EMCALClusterizer"), gBenchmark->GetCpuTime("EMCALClusterizer")/nEvents ); - } - + printf("Exec took %f seconds for Clusterizing", + gBenchmark->GetCpuTime("EMCALClusterizer")); + } } //____________________________________________________________________________ @@ -248,11 +283,12 @@ Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * emcRP, AliEMCALDigit ** for(iDigit = 0; iDigit < nDigits; iDigit++){ digit = maxAt[iDigit]; - Int_t relid[2] ; Float_t x = 0.; Float_t z = 0.; - fGeom->AbsToRelNumbering(digit->GetId(), relid) ; - fGeom->PosInAlice(relid, x, z) ; + // have to be tune for TRD1; May 31,06 + // Int_t relid[2] ; + // fGeom->AbsToRelNumbering(digit->GetId(), relid) ; // obsolete method + // fGeom->PosInAlice(relid, x, z) ; // obsolete method Float_t energy = maxAtEnergy[iDigit] ; @@ -316,19 +352,16 @@ void AliEMCALClusterizerv1::GetCalibrationParameters() // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB"); //Check if calibration is stored in data base - if(AliCDBManager::Instance()->IsDefaultStorageSet()){ - AliCDBEntry *entry = (AliCDBEntry*) AliCDBManager::Instance() - ->GetDefaultStorage() - ->Get("EMCAL/GainFactors_and_Pedestals/Calibration", - gAlice->GetRunNumber()); - if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject(); - } + + AliEMCALLoader *emcalLoader = + dynamic_cast(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")); + + fCalibData =emcalLoader->CalibData(); + if(!fCalibData) { //If calibration is not available use default parameters //Loader - AliEMCALLoader *emcalLoader = - dynamic_cast(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")); if ( !emcalLoader->Digitizer() ) emcalLoader->LoadDigitizer(); AliEMCALDigitizer * dig = dynamic_cast(emcalLoader->Digitizer()); @@ -345,7 +378,11 @@ void AliEMCALClusterizerv1::Init() // Attach the Clusterizer task to the list of EMCAL tasks AliRunLoader *rl = AliRunLoader::GetRunLoader(); - fGeom = dynamic_cast(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry(); + if (rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL")) + fGeom = dynamic_cast(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry(); + else + fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaulGeometryName()); + fGeom->GetTransformationForSM(); // Global <-> Local AliInfo(Form("geom 0x%x",fGeom)); @@ -363,11 +400,11 @@ void AliEMCALClusterizerv1::InitParameters() fNTowerInGroup = 36; //Produces maximum of 80 pseudoclusters per event - fECAClusteringThreshold = 0.5; // value obtained from Alexei + fECAClusteringThreshold = 0.1; // value obtained from Aleksei fECALocMaxCut = 0.03; // ?? - fECAW0 = 4.5 ; - fTimeGate = 1.e-8 ; + fECAW0 = 4.5; + fTimeCut = 300e-9 ; // 300 ns time cut (to be tuned) fToUnfold = kFALSE ; fRecPointsInRun = 0 ; fMinECut = 0.01; // have to be tune @@ -378,30 +415,31 @@ void AliEMCALClusterizerv1::InitParameters() //____________________________________________________________________________ Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2) const { - // Gives the neighbourness of two digits = 0 are not neighbour but continue searching + // Gives the neighbourness of two digits = 0 are not neighbour ; continue searching // = 1 are neighbour - // = 2 is in different SM, quit from searching + // = 2 is in different SM; continue searching // 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 static Int_t rv; - static Int_t nSupMod1=0, nTower1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0; - static Int_t nSupMod2=0, nTower2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0; + static Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0; + static Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0; static Int_t rowdiff, coldiff; rv = 0 ; - fGeom->GetCellIndex(d1->GetId(), nSupMod1,nTower1,nIphi1,nIeta1); - fGeom->GetCellIndex(d2->GetId(), nSupMod2,nTower2,nIphi2,nIeta2); + fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1); + fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2); if(nSupMod1 != nSupMod2) return 2; // different SM - fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nTower1,nIphi1,nIeta1, iphi1,ieta1); - fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nTower2,nIphi2,nIeta2, iphi2,ieta2); + fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1); + fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2); rowdiff = TMath::Abs(iphi1 - iphi2); coldiff = TMath::Abs(ieta1 - ieta2) ; - if (( coldiff <= 1 ) && ( rowdiff <= 1 )) rv = 1; // neighbours with at least commom vertex + // if (( coldiff <= 1 ) && ( rowdiff <= 1 )) rv = 1; // neighbours with at least commom vertex + if(coldiff + rowdiff <= 1) rv = 1; // neighbours with at least commom side; Nov 6,2006 if (gDebug == 2 && rv==1) printf("AreNeighbours: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d \n", @@ -429,21 +467,21 @@ Int_t AliEMCALClusterizerv1::AreInGroup(AliEMCALDigit * d1, AliEMCALDigit * d2) //find another in the list that will match? How about my TowerGroup search? static Int_t rv; - static Int_t nSupMod1=0, nTower1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0; - static Int_t nSupMod2=0, nTower2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0; + static Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0; + static Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0; static Int_t towerGroup1 = -1, towerGroup2 = -1; rv = 0 ; - fGeom->GetCellIndex(d1->GetId(), nSupMod1,nTower1,nIphi1,nIeta1); - fGeom->GetCellIndex(d2->GetId(), nSupMod2,nTower2,nIphi2,nIeta2); + fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1); + fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2); if(nSupMod1 != nSupMod2) return 2; // different SM - static Int_t nTowerInSM = fGeom->GetNCellsInSupMod()/fGeom->GetNCellsInTower(); + static Int_t nTowerInSM = fGeom->GetNCellsInSupMod()/fGeom->GetNCellsInModule(); //figure out which tower grouping each digit belongs to for(int it = 0; it < nTowerInSM/fNTowerInGroup; it++) { - if(nTower1 <= nTowerInSM - it*fNTowerInGroup) towerGroup1 = it; - if(nTower2 <= nTowerInSM - it*fNTowerInGroup) towerGroup2 = it; + if(nModule1 <= nTowerInSM - it*fNTowerInGroup) towerGroup1 = it; + if(nModule2 <= nTowerInSM - it*fNTowerInGroup) towerGroup2 = it; } if(towerGroup1 != towerGroup2) return 3; //different Towergroup @@ -457,16 +495,6 @@ Int_t AliEMCALClusterizerv1::AreInGroup(AliEMCALDigit * d1, AliEMCALDigit * d2) return rv ; } - -//____________________________________________________________________________ -void AliEMCALClusterizerv1::Unload() -{ - // Unloads the Digits and RecPoints - AliEMCALLoader *emcalLoader = dynamic_cast(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")); - - emcalLoader->UnloadDigits() ; - emcalLoader->UnloadRecPoints() ; -} //____________________________________________________________________________ void AliEMCALClusterizerv1::WriteRecPoints() @@ -485,11 +513,16 @@ void AliEMCALClusterizerv1::WriteRecPoints() emcalLoader->MakeRecPointsContainer(); treeR = emcalLoader->TreeR(); } + else if (treeR->GetEntries() > 0) { + Warning("WriteRecPoints","RecPoints already exist in output file. New Recpoitns will not be visible."); + } Int_t index ; //Evaluate position, dispersion and other RecPoint properties for EC section - for(index = 0; index < aECARecPoints->GetEntries(); index++) - (dynamic_cast(aECARecPoints->At(index)))->EvalAll(fECAW0,digits) ; + for(index = 0; index < aECARecPoints->GetEntries(); index++) { + if (dynamic_cast(aECARecPoints->At(index))->GetClusterType() != AliESDCaloCluster::kPseudoCluster) + dynamic_cast(aECARecPoints->At(index))->EvalAll(fECAW0,digits) ; + } aECARecPoints->Sort() ; @@ -511,7 +544,6 @@ void AliEMCALClusterizerv1::WriteRecPoints() treeR->Fill() ; emcalLoader->WriteRecPoints("OVERWRITE"); - } //____________________________________________________________________________ @@ -519,168 +551,122 @@ void AliEMCALClusterizerv1::MakeClusters(char* option) { // Steering method to construct the clusters stored in a list of Reconstructed Points // A cluster is defined as a list of neighbour digits + // Mar 03, 2007 by PAI - AliEMCALLoader *emcalLoader = dynamic_cast(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")); + if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader"); + AliEMCALLoader *emcalLoader = dynamic_cast(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")); TObjArray * aECARecPoints = emcalLoader->RecPoints() ; - - if (fGeom==0) - AliFatal("Did not get geometry from EMCALLoader"); - aECARecPoints->Clear(); - //Start with pseudoclusters, if option - if(strstr(option,"pseudo")) { - TClonesArray * digits = emcalLoader->Digits() ; - TClonesArray * digitsC = dynamic_cast(digits->Clone()); + TClonesArray *digits = emcalLoader->Digits(); - TIter nextdigit(digitsC) ; - AliEMCALDigit * digit; + // Set up TObjArray with pointers to digits to work on + TObjArray *digitsC = new TObjArray(); + TIter nextdigit(digits); + AliEMCALDigit *digit; + while ( (digit = dynamic_cast(nextdigit())) ) { + digitsC->AddLast(digit); + } - AliEMCALRecPoint * recPoint = 0 ; - int ineb=0; - nextdigit.Reset(); + //Start with pseudoclusters, if option + if(strstr(option,"pseudo")) { + // + // New algorithm : will be created one pseudo cluster per module + // + AliDebug(1,Form("Pseudo clustering #digits : %i ",digits->GetEntries())); - // PseudoClusterization starts - while ( (digit = dynamic_cast(nextdigit())) ) { // scan over the list of digitsC - recPoint = 0 ; - TArrayI clusterECAdigitslist(1000); // what is this + AliEMCALRecPoint *recPoints[12]; // max size is 12 : see fGeom->GetNumberOfSuperModules(); + for(int i=0; i<12; i++) recPoints[i] = 0; + TIter nextdigitC(digitsC) ; + // PseudoClusterization starts + int nSM = 0; // # of SM + while ( (digit = dynamic_cast(nextdigitC())) ) { // scan over the list of digitsC if(fGeom->CheckAbsCellId(digit->GetId()) ) { //Is this an EMCAL digit? Just maing sure... - Int_t iDigitInECACluster = 0; // start a new Tower RecPoint - if(fNumberOfECAClusters >= aECARecPoints->GetSize()) aECARecPoints->Expand(2*fNumberOfECAClusters+1) ; - AliEMCALRecPoint * rp = new AliEMCALRecPoint("") ; - aECARecPoints->AddAt(rp, fNumberOfECAClusters) ; - recPoint = dynamic_cast(aECARecPoints->At(fNumberOfECAClusters)) ; - fNumberOfECAClusters++ ; - - recPoint->SetClusterType(AliEMCALRecPoint::kPseudoCluster); - - recPoint->AddDigit(*digit, digit->GetAmp()) ; - clusterECAdigitslist[iDigitInECACluster] = digit->GetIndexInList() ; - iDigitInECACluster++ ; - digitsC->Remove(digit) ; - AliDebug(1,Form("MakePseudoClusters: OK id = %d, adc = %f \n", digit->GetId(), digit->GetAmp())); - nextdigit.Reset(); // will start from beggining - - AliEMCALDigit * digitN = 0; // digi neighbor - Int_t index = 0 ; - - // Find the neighbours - while (index < iDigitInECACluster){ // scan over digits already in cluster - digit = (AliEMCALDigit*)digits->At(clusterECAdigitslist[index]); - index++ ; - while ( (digitN = (AliEMCALDigit *)nextdigit())) { // scan over the reduced list of digits - ineb = AreInGroup(digit, digitN); // call (digit,digitN) in THAT oder !!!!! - switch (ineb ) { - case 0 : // not a neighbour - break ; - case 1 : // are neighbours - recPoint->AddDigit(*digitN, digitN->GetAmp() ) ; - clusterECAdigitslist[iDigitInECACluster] = digitN->GetIndexInList() ; - iDigitInECACluster++ ; - digitsC->Remove(digitN) ; - break ; - case 2 : // different SM - break ; - case 3 : // different Tower group - break ; - } // switch - } // scan over the reduced list of digits - } // scan over digits already in cluster - nextdigit.Reset() ; // will start from beggining + nSM = fGeom->GetSuperModuleNumber(digit->GetId()); + if(recPoints[nSM] == 0) { + recPoints[nSM] = new AliEMCALRecPoint(Form("PC%2.2i", nSM)); + recPoints[nSM]->SetClusterType(AliESDCaloCluster::kPseudoCluster); + } + recPoints[nSM]->AddDigit(*digit, Calibrate(digit->GetAmp(), digit->GetId())); } } - if(recPoint) cout << "cl.e " << recPoint->GetEnergy() << endl; - delete digitsC ; + fNumberOfECAClusters = 0; + for(int i=0; iGetNumberOfSuperModules(); i++) { // put non empty rec.points to container + if(recPoints[i]) aECARecPoints->AddAt(recPoints[i], fNumberOfECAClusters++); + } + AliDebug(1,Form(" Number of PC %d ", fNumberOfECAClusters)); } - //Now do real clusters - TClonesArray * digits = emcalLoader->Digits() ; - TClonesArray * digitsC = dynamic_cast(digits->Clone()); - - TIter nextdigit(digitsC) ; - AliEMCALDigit * digit; + // + // Now do real clusters + // - AliEMCALRecPoint * recPoint = 0 ; - int ineb=0; - nextdigit.Reset(); + double e = 0.0, ehs = 0.0; + TIter nextdigitC(digitsC); - double e=0.0, ehs = 0.0; - while ( (digit = dynamic_cast(nextdigit())) ) { // clean up digits + while ( (digit = dynamic_cast(nextdigitC())) ) { // clean up digits e = Calibrate(digit->GetAmp(), digit->GetId()); AliEMCALHistoUtilities::FillH1(fHists, 10, digit->GetAmp()); AliEMCALHistoUtilities::FillH1(fHists, 11, e); - if(e < fMinECut ) digitsC->Remove(digit); - else ehs += e; - } - cout << " Number of digits " << digits->GetEntries() << " -> (e>" <GetEntries()<< " ehs "<(nextdigit())) ) { // scan over the list of digitsC - recPoint = 0 ; - TArrayI clusterECAdigitslist(1000); // what is this + if ( e < fMinECut || digit->GetTimeR() > fTimeCut ) + digitsC->Remove(digit); + else + ehs += e; + } + AliDebug(1,Form("MakeClusters: Number of digits %d -> (e %f), ehs %d\n", + digits->GetEntries(),fMinECut,ehs)); + + nextdigitC.Reset(); + + while ( (digit = dynamic_cast(nextdigitC())) ) { // scan over the list of digitsC + TArrayI clusterECAdigitslist(digits->GetEntries()); if(fGeom->CheckAbsCellId(digit->GetId()) && (Calibrate(digit->GetAmp(), digit->GetId()) > fECAClusteringThreshold ) ){ - Int_t iDigitInECACluster = 0; // start a new Tower RecPoint + // start a new Tower RecPoint if(fNumberOfECAClusters >= aECARecPoints->GetSize()) aECARecPoints->Expand(2*fNumberOfECAClusters+1) ; - AliEMCALRecPoint * rp = new AliEMCALRecPoint("") ; - aECARecPoints->AddAt(rp, fNumberOfECAClusters) ; + AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("") ; + aECARecPoints->AddAt(recPoint, fNumberOfECAClusters) ; recPoint = dynamic_cast(aECARecPoints->At(fNumberOfECAClusters)) ; fNumberOfECAClusters++ ; - recPoint->SetClusterType(AliEMCALRecPoint::kClusterv1); + recPoint->SetClusterType(AliESDCaloCluster::kClusterv1); recPoint->AddDigit(*digit, Calibrate(digit->GetAmp(), digit->GetId())) ; - clusterECAdigitslist[iDigitInECACluster] = digit->GetIndexInList() ; - iDigitInECACluster++ ; + TObjArray clusterDigits; + clusterDigits.AddLast(digit); digitsC->Remove(digit) ; - AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , thre = %f \n", digit->GetId(), + + AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , cell.th. = %f \n", digit->GetId(), Calibrate(digit->GetAmp(),digit->GetId()), fECAClusteringThreshold)); - nextdigit.Reset(); // will start from beggining - AliEMCALDigit * digitN = 0; // digi neighbor - Int_t index = 0 ; - - // Find the neighbours - while (index < iDigitInECACluster){ // scan over digits already in cluster - digit = (AliEMCALDigit*)digits->At(clusterECAdigitslist[index]); - index++ ; - while ( (digitN = (AliEMCALDigit *)nextdigit())) { // scan over the reduced list of digits - ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!! - switch (ineb ) { - case 0 : // not a neighbour - break ; - case 1 : // are neighbours - recPoint->AddDigit(*digitN, Calibrate(digitN->GetAmp(),digitN->GetId()) ) ; - clusterECAdigitslist[iDigitInECACluster] = digitN->GetIndexInList() ; - iDigitInECACluster++ ; - digitsC->Remove(digitN) ; - break ; - case 2 : // different SM - break ; - } // switch - } // scan over the reduced list of digits + // Grow cluster by finding neighbours + TIter nextClusterDigit(&clusterDigits); + while ( (digit = dynamic_cast(nextClusterDigit())) ) { // scan over digits in cluster + TIter nextdigitN(digitsC); + AliEMCALDigit *digitN = 0; // digi neighbor + while ( (digitN = (AliEMCALDigit *)nextdigitN()) ) { // scan over all digits to look for neighbours + if (AreNeighbours(digit, digitN)==1) { // call (digit,digitN) in THAT oder !!!!! + recPoint->AddDigit(*digitN, Calibrate(digitN->GetAmp(),digitN->GetId()) ) ; + clusterDigits.AddLast(digitN) ; + digitsC->Remove(digitN) ; + } // if(ineb==1) + } // scan over digits } // scan over digits already in cluster - nextdigit.Reset() ; // will start from beggining - } - } // while digit - if(recPoint) cout << "cl.e " << recPoint->GetEnergy() << endl; + if(recPoint) + AliDebug(2,Form("MakeClusters: %d digitd, energy %f \n", clusterDigits.GetEntries(), recPoint->GetEnergy())); + } // If seed found + } // while digit + delete digitsC ; AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,digits->GetEntriesFast())); } -//____________________________________________________________________________ void AliEMCALClusterizerv1::MakeUnfolding() const { Fatal("AliEMCALClusterizerv1::MakeUnfolding", "--> Unfolding not implemented") ; - } //____________________________________________________________________________ @@ -760,24 +746,30 @@ void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option) AliEMCALLoader *emcalLoader = dynamic_cast(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")); TObjArray * aECARecPoints = emcalLoader->RecPoints() ; - printf("PrintRecPoints: Clusterization result:") ; + if(strstr(option,"deb")) { + printf("PrintRecPoints: Clusterization result:") ; - printf("event # %d\n", emcalLoader->GetRunLoader()->GetEventNumber() ) ; - printf(" Found %d ECA Rec Points\n ", + printf("event # %d\n", emcalLoader->GetRunLoader()->GetEventNumber() ) ; + printf(" Found %d ECA Rec Points\n ", aECARecPoints->GetEntriesFast()) ; + } fRecPointsInRun += aECARecPoints->GetEntriesFast() ; if(strstr(option,"all")) { - Int_t index =0; - printf("\n-----------------------------------------------------------------------\n") ; - printf("Clusters in ECAL section\n") ; - printf("Index Ene(GeV) Multi Module phi r theta X Y Z Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ; + if(strstr(option,"deb")) { + printf("\n-----------------------------------------------------------------------\n") ; + printf("Clusters in ECAL section\n") ; + printf("Index Ene(GeV) Multi Module GX GY GZ lX lY lZ Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ; + } + Int_t index =0; Float_t maxE=0; Float_t maxL1=0; Float_t maxL2=0; Float_t maxDis=0; + AliEMCALHistoUtilities::FillH1(fHists, 12, double(aECARecPoints->GetEntries())); + for (index = 0 ; index < aECARecPoints->GetEntries() ; index++) { AliEMCALRecPoint * rp = dynamic_cast(aECARecPoints->At(index)) ; TVector3 globalpos; @@ -789,6 +781,7 @@ void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option) Int_t * primaries; Int_t nprimaries; primaries = rp->GetPrimaries(nprimaries); + if(strstr(option,"deb")) printf("\n%6d %8.4f %3d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ", rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(), @@ -806,9 +799,11 @@ void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option) fPointDis->Fill(rp->GetDispersion()); fPointMult->Fill(rp->GetMultiplicity()); ///////////// - for (Int_t iprimary=0; iprimaryFill(maxE); @@ -816,27 +811,30 @@ void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option) fMaxL2->Fill(maxL2); fMaxDis->Fill(maxDis); - + if(strstr(option,"deb")) printf("\n-----------------------------------------------------------------------\n"); } } TList* AliEMCALClusterizerv1::BookHists() { + //set up histograms for monitoring clusterizer performance + gROOT->cd(); - fPointE = new TH1F("pointE","point energy", 2000, 0.0, 150.); - fPointL1 = new TH1F("pointL1","point L1", 1000, 0.0, 3.); - fPointL2 = new TH1F("pointL2","point L2", 1000, 0.0, 3.); - fPointDis = new TH1F("pointDis","point Dis", 1000, 0.0, 3.); - fPointMult = new TH1F("pointMult","point Mult", 100, 0.0, 100.); - fDigitAmp = new TH1F("digitAmp","Digit Amplitude", 2000, 0.0, 5000.); - fMaxE = new TH1F("maxE","Max point energy", 2000, 0.0, 150.); - fMaxL1 = new TH1F("maxL1","Max point L1", 1000, 0.0, 3.); - fMaxL2 = new TH1F("maxL2","Max point L2", 1000, 0.0, 3.); - fMaxDis = new TH1F("maxDis","Max point Dis", 1000, 0.0, 3.); // 9 + fPointE = new TH1F("00_pointE","point energy", 2000, 0.0, 150.); + fPointL1 = new TH1F("01_pointL1","point L1", 1000, 0.0, 3.); + fPointL2 = new TH1F("02_pointL2","point L2", 1000, 0.0, 3.); + fPointDis = new TH1F("03_pointDisp","point dispersion", 1000, 0.0, 10.); + fPointMult = new TH1F("04_pointMult","#cell in point(cluster)", 101, -0.5, 100.5); + fDigitAmp = new TH1F("05_digitAmp","Digit Amplitude", 2000, 0.0, 5000.); + fMaxE = new TH1F("06_maxE","Max point energy", 2000, 0.0, 150.); + fMaxL1 = new TH1F("07_maxL1","Largest (first) of eigenvalue of covariance matrix", 1000, 0.0, 3.); + fMaxL2 = new TH1F("08_maxL2","Smalest (second) of eigenvalue of covariace matrix", 1000, 0.0, 3.); + fMaxDis = new TH1F("09_maxDis","Point dispersion", 1000, 0.0, 10.); // 9 // - new TH1F("adcOfDigits","adc of digits(threshold control)", 1001, -0.5, 1000.5); // 10 - new TH1F("energyOfDigits","energy of digits(threshold control)", 1000, 0.0, 1.); // 11 + new TH1F("10_adcOfDigits","adc of digits(threshold control)", 1001, -0.5, 1000.5); // 10 + new TH1F("11_energyOfDigits","energy of digits(threshold control)", 1000, 0.0, 1.); // 11 + new TH1F("12_numberOfPoints","number of points(clusters)", 101, -0.5, 100.5); // 12 return AliEMCALHistoUtilities::MoveHistsToList("EmcalClusterizerv1ControlHists", kFALSE); } @@ -846,8 +844,34 @@ void AliEMCALClusterizerv1::SaveHists(const char *fn) AliEMCALHistoUtilities::SaveListOfHists(fHists, fn, kTRUE); } +void AliEMCALClusterizerv1::PrintRecoInfo() +{ + printf(" AliEMCALClusterizerv1::PrintRecoInfo() : version %s \n", Version() ); + TH1F *h = (TH1F*)fHists->At(12); + if(h) { + printf(" ## Multiplicity of RecPoints ## \n"); + for(int i=1; i<=h->GetNbinsX(); i++) { + int nbin = int((*h)[i]); + int mult = int(h->GetBinCenter(i)); + if(nbin > 0) printf(" %i : %5.5i %6.3f %% \n", mult, nbin, 100.*nbin/h->GetEntries()); + } + } +} + +void AliEMCALClusterizerv1::DrawLambdasHists() +{ + if(fMaxL1) { + fMaxL1->Draw(); + if(fMaxL2) fMaxL2->Draw("same"); + if(fMaxDis) { + fMaxDis->Draw("same"); + } + } +} + void AliEMCALClusterizerv1::Browse(TBrowser* b) { if(fHists) b->Add(fHists); + if(fGeom) b->Add(fGeom); TTask::Browse(b); }