X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliEMCALClusterizerv1.cxx;h=00c43d9cce97731fc79d8a89292f138d1da99c11;hb=12fb44337a22da7b068695e74a8e8e529a141e38;hp=63a187819723ced2bbf62042431f0783e34a753a;hpb=70a931984a8b03dd027bab16bbe2ae14989707f8;p=u%2Fmrichter%2FAliRoot.git diff --git a/EMCAL/AliEMCALClusterizerv1.cxx b/EMCAL/AliEMCALClusterizerv1.cxx index 63a18781972..00c43d9cce9 100644 --- a/EMCAL/AliEMCALClusterizerv1.cxx +++ b/EMCAL/AliEMCALClusterizerv1.cxx @@ -15,9 +15,10 @@ /* $Id$ */ -//*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute) +//-- 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,122 +46,160 @@ // // time - print benchmarking results // --- ROOT system --- - -#include "TROOT.h" -#include "TFile.h" -#include "TFolder.h" -#include "TMath.h" -#include "TMinuit.h" -#include "TTree.h" -#include "TSystem.h" -#include "TBenchmark.h" +#include + +class TROOT; +#include +#include +class TFolder; +#include +#include +#include +class TSystem; +#include +#include +#include // --- Standard library --- // --- AliRoot header files --- -#include "AliEMCALGetter.h" +#include "AliRunLoader.h" +#include "AliRun.h" +#include "AliESD.h" #include "AliEMCALClusterizerv1.h" #include "AliEMCALRecPoint.h" #include "AliEMCALDigit.h" #include "AliEMCALDigitizer.h" #include "AliEMCAL.h" #include "AliEMCALGeometry.h" +#include "AliEMCALHistoUtilities.h" +#include "AliEMCALRecParam.h" +#include "AliEMCALReconstructor.h" +#include "AliCDBManager.h" + +class AliCDBStorage; +#include "AliCDBEntry.h" ClassImp(AliEMCALClusterizerv1) - -//____________________________________________________________________________ - AliEMCALClusterizerv1::AliEMCALClusterizerv1() : AliEMCALClusterizer() -{ - // default ctor (to be used mainly by Streamer) - - InitParameters() ; - fDefaultInit = kTRUE ; -} //____________________________________________________________________________ -AliEMCALClusterizerv1::AliEMCALClusterizerv1(const TString alirunFileName, const TString eventFolderName) -:AliEMCALClusterizer(alirunFileName, eventFolderName) +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(kFALSE), + fToUnfold(kFALSE), + fNumberOfECAClusters(0),fNTowerInGroup(0),fCalibData(0), + fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.), + fECAW0(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() { // dtor - } //____________________________________________________________________________ -const TString AliEMCALClusterizerv1::BranchName() const -{ - return GetName(); +Float_t AliEMCALClusterizerv1::Calibrate(Int_t amp, Int_t AbsId) +{ + + // Convert digitized amplitude into energy. + // Calibration parameters are taken from calibration data base for raw data, + // or from digitizer parameters for simulated data. -} + if(fCalibData){ + + if (fGeom==0) + AliFatal("Did not get geometry from EMCALLoader") ; + + Int_t iSupMod = -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, nModule, nIphi, nIeta) ; + if(!bCell) { + fGeom->PrintGeometry(); + Error("Calibrate()"," Wrong cell id number : %i", AbsId); + assert(0); + } -//____________________________________________________________________________ -Float_t AliEMCALClusterizerv1::Calibrate(Int_t amp) const -{ - //To be replased later by the method, reading individual parameters from the database - return -fADCpedestalECA + amp * fADCchannelECA ; + fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta); + + fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi); + fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi); + + return -fADCpedestalECA + amp * fADCchannelECA ; + + } + else //Return energy with default parameters if calibration is not available + return -fADCpedestalECA + amp * fADCchannelECA ; + } //____________________________________________________________________________ -void AliEMCALClusterizerv1::Exec(Option_t * option) +void AliEMCALClusterizerv1::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 the current event + // in AliEMCALLoader if(strstr(option,"tim")) gBenchmark->Start("EMCALClusterizer"); if(strstr(option,"print")) Print("") ; + + //Get calibration parameters from file or digitizer default values. + GetCalibrationParameters() ; - AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle()) ; - - if (fLastEvent == -1) - fLastEvent = gime->MaxEvent() - 1; - - Int_t nEvents = fLastEvent - fFirstEvent + 1; - Int_t ievent ; - for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) { - gime->Event(ievent,"D") ; + fNumberOfECAClusters = 0; - GetCalibrationParameters() ; + if(strstr(option,"pseudo")) + MakeClusters("pseudo") ; //both types + else + MakeClusters("") ; //only the real clusters - fNumberOfECAClusters = 0; - - MakeClusters() ; + if(fToUnfold) + MakeUnfolding() ; - if(fToUnfold) - MakeUnfolding() ; + Int_t index ; - WriteRecPoints() ; + //Evaluate position, dispersion and other RecPoint properties for EC section + for(index = 0; index < fRecPoints->GetEntries(); index++) { + if (dynamic_cast(fRecPoints->At(index))->GetClusterType() != AliESDCaloCluster::kEMCALPseudoCluster) + dynamic_cast(fRecPoints->At(index))->EvalAll(fECAW0,fDigitsArr) ; + } - if(strstr(option,"deb")) - PrintRecPoints(option) ; + fRecPoints->Sort() ; - //increment the total number of recpoints per run - fRecPointsInRun += gime->ECARecPoints()->GetEntriesFast() ; + for(index = 0; index < fRecPoints->GetEntries(); index++) { + (dynamic_cast(fRecPoints->At(index)))->SetIndexInList(index) ; + (dynamic_cast(fRecPoints->At(index)))->Print(); } + + fTreeR->Fill(); - Unload(); + if(strstr(option,"deb") || strstr(option,"all")) + PrintRecPoints(option) ; + + AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->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")); + } } //____________________________________________________________________________ @@ -171,16 +210,13 @@ Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * emcRP, AliEMCALDigit ** // The initial values for fitting procedure are set equal to the positions of local maxima. // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers - AliEMCALGetter * gime = AliEMCALGetter::Instance() ; - TClonesArray * digits = gime->Digits() ; - gMinuit->mncler(); // Reset Minuit's list of paramters gMinuit->SetPrintLevel(-1) ; // No Printout gMinuit->SetFCN(AliEMCALClusterizerv1::UnfoldingChiSquare) ; // To set the address of the minimization function TList * toMinuit = new TList(); toMinuit->AddAt(emcRP,0) ; - toMinuit->AddAt(digits,1) ; + toMinuit->AddAt(fDigitsArr,1) ; gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare @@ -193,16 +229,15 @@ Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * emcRP, AliEMCALDigit ** Int_t iDigit ; - AliEMCALGeometry * geom = gime->EMCALGeometry() ; - for(iDigit = 0; iDigit < nDigits; iDigit++){ digit = maxAt[iDigit]; - Int_t relid[2] ; Float_t x = 0.; Float_t z = 0.; - geom->AbsToRelNumbering(digit->GetId(), relid) ; - geom->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] ; @@ -257,15 +292,38 @@ Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * emcRP, AliEMCALDigit ** //____________________________________________________________________________ void AliEMCALClusterizerv1::GetCalibrationParameters() { - // Gets the parameters for the calibration from the digitizer - AliEMCALGetter * gime = AliEMCALGetter::Instance() ; - - if ( !gime->Digitizer() ) - gime->LoadDigitizer(); - AliEMCALDigitizer * dig = gime->Digitizer(); - - fADCchannelECA = dig->GetECAchannel() ; - fADCpedestalECA = dig->GetECApedestal(); + // Set calibration parameters: + // if calibration database exists, they are read from database, + // otherwise, they are taken from digitizer. + // + // It is a user responsilibity to open CDB before reconstruction, + // for example: + // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB"); + + //Check if calibration is stored in data base + + if(!fCalibData && (AliCDBManager::Instance()->IsDefaultStorageSet())) + { + AliCDBEntry *entry = (AliCDBEntry*) + AliCDBManager::Instance()->Get("EMCAL/Calib/Data"); + if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject(); + } + + if(!fCalibData) + AliFatal("Calibration parameters not found in CDB!"); + + // Please fix it!! Or better just remove it... +// if(!fCalibData) +// { +// //If calibration is not available use default parameters +// //Loader +// if ( !emcalLoader->Digitizer() ) +// emcalLoader->LoadDigitizer(); +// AliEMCALDigitizer * dig = dynamic_cast(emcalLoader->Digitizer()); + +// fADCchannelECA = dig->GetECAchannel() ; +// fADCpedestalECA = dig->GetECApedestal(); +// } } //____________________________________________________________________________ @@ -274,16 +332,18 @@ void AliEMCALClusterizerv1::Init() // Make all memory allocations which can not be done in default constructor. // Attach the Clusterizer task to the list of EMCAL tasks - AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle(), fEventFolderName.Data()); + AliRunLoader *rl = AliRunLoader::GetRunLoader(); + if (rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL")) + fGeom = dynamic_cast(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry(); + else + fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaulGeometryName()); - AliEMCALGeometry * geom = gime->EMCALGeometry() ; + AliDebug(1,Form("geom 0x%x",fGeom)); - fNTowers = geom->GetNZ() * geom->GetNPhi() ; if(!gMinuit) gMinuit = new TMinuit(100) ; - if ( !gime->Clusterizer() ) - gime->PostClusterizer(this); + fHists = BookHists(); } //____________________________________________________________________________ @@ -291,212 +351,231 @@ void AliEMCALClusterizerv1::InitParameters() { // Initializes the parameters for the Clusterizer fNumberOfECAClusters = 0; - fECAClusteringThreshold = 0.0045; // must be adjusted according to the noise leve set by digitizer - fECALocMaxCut = 0.03 ; - fECAW0 = 4.5 ; - fTimeGate = 1.e-8 ; + + fNTowerInGroup = 36; //Produces maximum of 80 pseudoclusters per event + + fECALocMaxCut = 0.03; // ?? + + fTimeCut = 300e-9 ; // 300 ns time cut (to be tuned) fToUnfold = kFALSE ; - fRecPointsInRun = 0 ; + + fCalibData = 0 ; + + const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam(); + if(!recParam) { + AliFatal("Reconstruction parameters for EMCAL not set!"); + } + else { + fECAClusteringThreshold = recParam->GetClusteringThreshold(); + fECAW0 = recParam->GetW0(); + fMinECut = recParam->GetMinECut(); + AliDebug(1,Form("Reconstruction parameters: fECAClusteringThreshold=%.3f, fECAW=%.3f, fMinECut=%.3f", + fECAClusteringThreshold,fECAW0,fMinECut)); + } + } //____________________________________________________________________________ -Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2)const +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 are not neighbour but do not continue 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 - AliEMCALGeometry * geom = AliEMCALGetter::Instance()->EMCALGeometry() ; + static Int_t rv; + 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 ; - Int_t rv = 0 ; + fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1); + fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2); + if(nSupMod1 != nSupMod2) return 2; // different SM - Int_t relid1[2] ; - geom->AbsToRelNumbering(d1->GetId(), relid1) ; - - Int_t relid2[2] ; - geom->AbsToRelNumbering(d2->GetId(), relid2) ; - + fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1); + fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2); - Int_t rowdiff = TMath::Abs( relid1[0] - relid2[0] ) ; - Int_t coldiff = TMath::Abs( relid1[1] - relid2[1] ) ; + rowdiff = TMath::Abs(iphi1 - iphi2); + coldiff = TMath::Abs(ieta1 - ieta2) ; - if (( coldiff <= 1 ) && ( rowdiff <= 1 )){ - if(TMath::Abs(d1->GetTime() - d2->GetTime() ) < fTimeGate) - rv = 1 ; - } - else { - if((relid2[0] > relid1[0]) && (relid2[1] > relid1[1]+1)) - rv = 2; // Difference in row numbers is too large to look further - } + // neighbours with at least commom side; May 11, 2007 + if ((coldiff==0 && abs(rowdiff)==1) || (rowdiff==0 && abs(coldiff)==1)) rv = 1; - if (gDebug == 2 ) - printf("AreNeighbours: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d ", - rv, d1->GetId(), relid1[0], relid1[1], - d2->GetId(), relid2[0], relid2[1]) ; + if (gDebug == 2 && rv==1) + printf("AreNeighbours: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d \n", + rv, d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2); return rv ; } //____________________________________________________________________________ -void AliEMCALClusterizerv1::Unload() +Int_t AliEMCALClusterizerv1::AreInGroup(AliEMCALDigit * d1, AliEMCALDigit * d2) const { - // Unloads the Digits and RecPoints - AliEMCALGetter * gime = AliEMCALGetter::Instance() ; - gime->EmcalLoader()->UnloadDigits() ; - gime->EmcalLoader()->UnloadRecPoints() ; -} - -//____________________________________________________________________________ -void AliEMCALClusterizerv1::WriteRecPoints() -{ - - // Creates new branches with given title - // fills and writes into TreeR. - - AliEMCALGetter *gime = AliEMCALGetter::Instance() ; - - TObjArray * aECARecPoints = gime->ECARecPoints() ; - - TClonesArray * digits = gime->Digits() ; - TTree * treeR = gime->TreeR(); ; - - 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) ; - - aECARecPoints->Sort() ; - - for(index = 0; index < aECARecPoints->GetEntries(); index++) - (dynamic_cast(aECARecPoints->At(index)))->SetIndexInList(index) ; - - aECARecPoints->Expand(aECARecPoints->GetEntriesFast()) ; - - Int_t bufferSize = 32000 ; - Int_t splitlevel = 0 ; + // Tells whether two digits fall within the same supermodule and + // tower grouping. The number of towers in a group is controlled by + // the parameter nTowersInGroup + // = 0 are not in same group but continue searching + // = 1 same group + // = 2 is in different SM, quit from searching + // = 3 different tower group, quit from searching + // + // 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 + + //JLK Question: does the quit from searching assume that the digits + //are ordered, so that once you are in a different SM, you'll not + //find another in the list that will match? How about my TowerGroup search? + + static Int_t rv; + 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,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->GetNCellsInModule(); + + //figure out which tower grouping each digit belongs to + for(int it = 0; it < nTowerInSM/fNTowerInGroup; it++) { + if(nModule1 <= nTowerInSM - it*fNTowerInGroup) towerGroup1 = it; + if(nModule2 <= nTowerInSM - it*fNTowerInGroup) towerGroup2 = it; + } + if(towerGroup1 != towerGroup2) return 3; //different Towergroup - //EC section branch - TBranch * branchECA = treeR->Branch("EMCALECARP","TObjArray",&aECARecPoints,bufferSize,splitlevel); - branchECA->SetTitle(BranchName()); + //same SM, same towergroup, we're happy + if(towerGroup1 == towerGroup2 && towerGroup2 >= 0) + rv = 1; - branchECA->Fill() ; + if (gDebug == 2 && rv==1) + printf("AreInGroup: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d \n", + rv, d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2); - gime->WriteRecPoints("OVERWRITE"); - gime->WriteClusterizer("OVERWRITE"); + return rv ; } - + //____________________________________________________________________________ -void AliEMCALClusterizerv1::MakeClusters() +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 - - AliEMCALGetter * gime = AliEMCALGetter::Instance() ; + // Mar 03, 2007 by PAI - AliEMCALGeometry * geom = gime->EMCALGeometry() ; + if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader"); - TObjArray * aECARecPoints = gime->ECARecPoints() ; + fRecPoints->Clear(); - aECARecPoints->Delete() ; - - TClonesArray * digits = gime->Digits() ; - - TIter next(digits) ; - AliEMCALDigit * digit ; - Int_t ndigECA=0 ; + // Set up TObjArray with pointers to digits to work on + TObjArray *digitsC = new TObjArray(); + TIter nextdigit(fDigitsArr); + AliEMCALDigit *digit; + while ( (digit = dynamic_cast(nextdigit())) ) { + digitsC->AddLast(digit); + } - // count the number of digits in ECA section - while ( (digit = dynamic_cast(next())) ) { // scan over the list of digits - if (geom->IsInECA(digit->GetId())) - ndigECA++ ; - else { - Error("MakeClusters", "id = %d is a wrong ID!", digit->GetId()) ; - abort() ; + //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 ",fDigitsArr->GetEntries())); + + 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... + nSM = fGeom->GetSuperModuleNumber(digit->GetId()); + if(recPoints[nSM] == 0) { + recPoints[nSM] = new AliEMCALRecPoint(Form("PC%2.2i", nSM)); + recPoints[nSM]->SetClusterType(AliESDCaloCluster::kEMCALPseudoCluster); + } + recPoints[nSM]->AddDigit(*digit, Calibrate(digit->GetAmp(), digit->GetId())); + } } - } - TClonesArray * digitsC = dynamic_cast(digits->Clone()) ; - // Clusterization starts - - TIter nextdigit(digitsC) ; - - while ( (digit = dynamic_cast(nextdigit())) ) { // scan over the list of digitsC - AliEMCALRecPoint * clu = 0 ; - - TArrayI clusterECAdigitslist(50); - - Bool_t inECA = kFALSE; - if( geom->IsInECA(digit->GetId()) ) { - inECA = kTRUE ; - } - if (gDebug == 2) { - if (inECA) - printf("MakeClusters: id = %d, ene = %f , thre = %f", - digit->GetId(),Calibrate(digit->GetAmp()), fECAClusteringThreshold) ; + fNumberOfECAClusters = 0; + for(int i=0; iGetNumberOfSuperModules(); i++) { // put non empty rec.points to container + if(recPoints[i]) fRecPoints->AddAt(recPoints[i], fNumberOfECAClusters++); } - if (inECA && (Calibrate(digit->GetAmp()) > fECAClusteringThreshold ) ){ - - Int_t iDigitInECACluster = 0; - // Find the seed - if( geom->IsInECA(digit->GetId()) ) { - // start a new Tower RecPoint - if(fNumberOfECAClusters >= aECARecPoints->GetSize()) - aECARecPoints->Expand(2*fNumberOfECAClusters+1) ; - AliEMCALRecPoint * rp = new AliEMCALRecPoint("") ; - aECARecPoints->AddAt(rp, fNumberOfECAClusters) ; - clu = dynamic_cast(aECARecPoints->At(fNumberOfECAClusters)) ; - fNumberOfECAClusters++ ; - clu->AddDigit(*digit, Calibrate(digit->GetAmp())) ; - clusterECAdigitslist[iDigitInECACluster] = digit->GetIndexInList() ; - iDigitInECACluster++ ; - digitsC->Remove(digit) ; - if (gDebug == 2 ) - printf("MakeClusters: OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp()), fECAClusteringThreshold) ; - - } - nextdigit.Reset() ; + AliDebug(1,Form(" Number of PC %d ", fNumberOfECAClusters)); + } + + // + // Now do real clusters + // + + double e = 0.0, ehs = 0.0; + TIter nextdigitC(digitsC); + + 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 || digit->GetTimeR() > fTimeCut ) + digitsC->Remove(digit); + else + ehs += e; + } + AliDebug(1,Form("MakeClusters: Number of digits %d -> (e %f), ehs %d\n", + fDigitsArr->GetEntries(),fMinECut,ehs)); + + nextdigitC.Reset(); + + while ( (digit = dynamic_cast(nextdigitC())) ) { // scan over the list of digitsC + TArrayI clusterECAdigitslist(fDigitsArr->GetEntries()); + + if(fGeom->CheckAbsCellId(digit->GetId()) && (Calibrate(digit->GetAmp(), digit->GetId()) > fECAClusteringThreshold ) ){ + // start a new Tower RecPoint + if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(2*fNumberOfECAClusters+1) ; + AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("") ; + fRecPoints->AddAt(recPoint, fNumberOfECAClusters) ; + recPoint = dynamic_cast(fRecPoints->At(fNumberOfECAClusters)) ; + fNumberOfECAClusters++ ; + + recPoint->SetClusterType(AliESDCaloCluster::kEMCALClusterv1); + + recPoint->AddDigit(*digit, Calibrate(digit->GetAmp(), digit->GetId())) ; + TObjArray clusterDigits; + clusterDigits.AddLast(digit); + digitsC->Remove(digit) ; + + AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , cell.th. = %f \n", digit->GetId(), + Calibrate(digit->GetAmp(),digit->GetId()), fECAClusteringThreshold)); - AliEMCALDigit * digitN ; - Int_t index = 0 ; - - // Do the Clustering - - 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 - Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!! - switch (ineb ) { - case 0 : // not a neighbour - break ; - case 1 : // are neighbours - clu->AddDigit(*digitN, Calibrate( digitN->GetAmp()) ) ; - clusterECAdigitslist[iDigitInECACluster] = digitN->GetIndexInList() ; - iDigitInECACluster++ ; - digitsC->Remove(digitN) ; - break ; - case 2 : // too far from each other - goto endofloop1; - } // switch - - } // while digitN - - endofloop1: ; - nextdigit.Reset() ; - } // loop over ECA cluster - } // energy theshold - } // while digit + // 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 + 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,fDigitsArr->GetEntriesFast())); } -//____________________________________________________________________________ void AliEMCALClusterizerv1::MakeUnfolding() const { Fatal("AliEMCALClusterizerv1::MakeUnfolding", "--> Unfolding not implemented") ; - } //____________________________________________________________________________ @@ -544,18 +623,12 @@ void AliEMCALClusterizerv1::Print(Option_t * /*option*/)const // Print parameters - TString taskName(GetName()) ; - taskName.ReplaceAll(Version(), "") ; + TString taskName(Version()) ; printf("--------------- "); printf(taskName.Data()) ; printf(" "); - printf(GetTitle()) ; - printf("-----------\n"); - printf("Clusterizing digits from the file: "); - printf(taskName.Data()); - printf("\n Branch: "); - printf(GetName()); + printf("Clusterizing digits: "); printf("\n ECA Local Maximum cut = %f", fECALocMaxCut); printf("\n ECA Logarithmic weight = %f", fECAW0); if(fToUnfold) @@ -573,26 +646,31 @@ void AliEMCALClusterizerv1::Print(Option_t * /*option*/)const void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option) { // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer - - TObjArray * aECARecPoints = AliEMCALGetter::Instance()->ECARecPoints() ; - printf("PrintRecPoints: Clusterization result:") ; + if(strstr(option,"deb")) { + printf("PrintRecPoints: Clusterization result:") ; - printf("event # %d\n", gAlice->GetEvNumber() ) ; - printf(" Found %d ECA Rec Points\n ", - aECARecPoints->GetEntriesFast()) ; + printf(" Found %d ECA Rec Points\n ", + fRecPoints->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") ; - - for (index = 0 ; index < aECARecPoints->GetEntries() ; index++) { - AliEMCALRecPoint * rp = dynamic_cast(aECARecPoints->At(index)) ; + 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(fRecPoints->GetEntries())); + + for (index = 0 ; index < fRecPoints->GetEntries() ; index++) { + AliEMCALRecPoint * rp = dynamic_cast(fRecPoints->At(index)) ; TVector3 globalpos; - rp->GetGlobalPosition(globalpos); + //rp->GetGlobalPosition(globalpos); TVector3 localpos; rp->GetLocalPosition(localpos); Float_t lambda[2]; @@ -600,14 +678,90 @@ 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(), rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ; - for (Int_t iprimary=0; iprimaryGetEnergy()>maxE){ + maxE=rp->GetEnergy(); + maxL1=lambda[0]; + maxL2=lambda[1]; + maxDis=rp->GetDispersion(); + } + fPointE->Fill(rp->GetEnergy()); + fPointL1->Fill(lambda[0]); + fPointL2->Fill(lambda[1]); + fPointDis->Fill(rp->GetDispersion()); + fPointMult->Fill(rp->GetMultiplicity()); + ///////////// + if(strstr(option,"deb")){ + for (Int_t iprimary=0; iprimaryFill(maxE); + fMaxL1->Fill(maxL1); + 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("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("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); +} + +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"); + } + } +}