X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliEMCALClusterizerv1.cxx;h=044b6b334dbfc97f051f8ec3878c520f20da8bb5;hb=b97c26da4b2d527d08a1b07186508293b98604c3;hp=2daf060bf58fbd264d85a5216c8321c3b5dc228c;hpb=ed611565254f91186cbbd2b53f22a523d3ed7b50;p=u%2Fmrichter%2FAliRoot.git diff --git a/EMCAL/AliEMCALClusterizerv1.cxx b/EMCAL/AliEMCALClusterizerv1.cxx index 2daf060bf58..044b6b334db 100644 --- a/EMCAL/AliEMCALClusterizerv1.cxx +++ b/EMCAL/AliEMCALClusterizerv1.cxx @@ -14,16 +14,7 @@ **************************************************************************/ /* $Id$ */ -/* $Log: - 1 October 2000. Yuri Kharlov: - AreNeighbours() - PPSD upper layer is considered if number of layers>1 - 18 October 2000. Yuri Kharlov: - AliEMCALClusterizerv1() - CPV clusterizing parameters added - MakeClusters() - After first PPSD digit remove EMC digits only once -*/ + //*-- 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) @@ -55,91 +46,127 @@ // --- 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" +class TROOT; +#include +#include +class TFolder; +#include +#include +#include +class TSystem; +#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" #include "AliEMCALDigit.h" #include "AliEMCALDigitizer.h" -#include "AliEMCALTowerRecPoint.h" #include "AliEMCAL.h" -#include "AliEMCALGetter.h" #include "AliEMCALGeometry.h" -#include "AliRun.h" +#include "AliEMCALHistoUtilities.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 ; + fDefaultInit = kTRUE ; + fGeom = AliEMCALGeometry::GetInstance(); + fGeom->GetTransformationForSM(); // Global <-> Local } //____________________________________________________________________________ -AliEMCALClusterizerv1::AliEMCALClusterizerv1(const char* headerFile, const char* name, const Bool_t toSplit) -:AliEMCALClusterizer(headerFile, name, toSplit) +AliEMCALClusterizerv1::AliEMCALClusterizerv1(const TString alirunFileName, const TString eventFolderName) +:AliEMCALClusterizer(alirunFileName, eventFolderName) { // ctor with the indication of the file where header Tree and digits Tree are stored InitParameters() ; Init() ; fDefaultInit = kFALSE ; - } //____________________________________________________________________________ AliEMCALClusterizerv1::~AliEMCALClusterizerv1() { // dtor - fSplitFile = 0 ; - } //____________________________________________________________________________ const TString AliEMCALClusterizerv1::BranchName() const { - TString branchName(GetName() ) ; - branchName.Remove(branchName.Index(Version())-1) ; - return branchName ; + return GetName(); } //____________________________________________________________________________ -Float_t AliEMCALClusterizerv1::Calibrate(Int_t amp, Int_t where) const +Float_t AliEMCALClusterizerv1::Calibrate(Int_t amp, Int_t AbsId) { - //To be replased later by the method, reading individual parameters from the database - // where = 0 == PRE ; where = 1 == ECAL ; where = 2 == HCAL - if ( where == 0 ) // calibrate as PRE section - return -fADCpedestalPRE + amp * fADCchannelPRE ; - else if (where == 1) //calibrate as EC section - return -fADCpedestalEC + amp * fADCchannelEC ; - else if (where == 2) //calibrate as HC section - return -fADCpedestalHC + amp * fADCchannelHC ; - else - Fatal("Calibrate", "Something went wrong!") ; - return -9999999. ; + + // 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){ + + //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 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); + + 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) { - // Steering method - - if( strcmp(GetName(), "")== 0 ) - Init() ; + // 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. if(strstr(option,"tim")) gBenchmark->Start("EMCALClusterizer"); @@ -147,56 +174,60 @@ void AliEMCALClusterizerv1::Exec(Option_t * option) if(strstr(option,"print")) Print("") ; - AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; - if(gime->BranchExists("RecPoints")) - return ; - Int_t nevents = gime->MaxEvent() ; - Int_t ievent ; + AliRunLoader *rl = AliRunLoader::GetRunLoader(); + AliEMCALLoader *emcalLoader = dynamic_cast(rl->GetDetectorLoader("EMCAL")); - for(ievent = 0; ievent < nevents; ievent++){ + //Get calibration parameters from file or digitizer default values. + GetCalibrationParameters() ; - gime->Event(ievent,"D") ; + if (fLastEvent == -1) + fLastEvent = rl->GetNumberOfEvents() - 1; + Int_t nEvents = fLastEvent - fFirstEvent + 1; - if(ievent == 0) - GetCalibrationParameters() ; + Int_t ievent ; + rl->LoadDigits("EMCAL"); + for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) { + rl->GetEvent(ievent); - fNumberOfPREClusters = fNumberOfECClusters = fNumberOfHCClusters = 0 ; - - MakeClusters() ; + fNumberOfECAClusters = 0; + + if(strstr(option,"pseudo")) + MakeClusters("pseudo") ; //both types + else + MakeClusters("") ; //only the real clusters if(fToUnfold) MakeUnfolding() ; - WriteRecPoints(ievent) ; + WriteRecPoints() ; if(strstr(option,"deb")) PrintRecPoints(option) ; - //increment the total number of digits per run - fRecPointsInRun += gime->PRERecPoints()->GetEntriesFast() ; - fRecPointsInRun += gime->ECALRecPoints()->GetEntriesFast() ; - fRecPointsInRun += gime->HCALRecPoints()->GetEntriesFast() ; - } + //increment the total number of recpoints per run + fRecPointsInRun += emcalLoader->RecPoints()->GetEntriesFast() ; + } + Unload(); + if(strstr(option,"tim")){ gBenchmark->Stop("EMCALClusterizer"); - Info("Exec", "took %f seconds for Clusterizing %f seconds per event", - gBenchmark->GetCpuTime("EMCALClusterizer"), gBenchmark->GetCpuTime("EMCALClusterizer")/nevents ) ; + printf("Exec took %f seconds for Clusterizing %f seconds per event", + gBenchmark->GetCpuTime("EMCALClusterizer"), gBenchmark->GetCpuTime("EMCALClusterizer")/nEvents ); } - + } //____________________________________________________________________________ -Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALTowerRecPoint * emcRP, AliEMCALDigit ** maxAt, Float_t * maxAtEnergy, +Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * emcRP, AliEMCALDigit ** maxAt, Float_t * maxAtEnergy, Int_t nPar, Float_t * fitparameters) const { // Calls TMinuit to fit the energy distribution of a cluster with several maxima // 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::GetInstance() ; - TClonesArray * digits = gime->Digits() ; - + AliEMCALLoader *emcalLoader = dynamic_cast(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")); + TClonesArray *digits = emcalLoader->Digits(); gMinuit->mncler(); // Reset Minuit's list of paramters gMinuit->SetPrintLevel(-1) ; // No Printout @@ -217,16 +248,15 @@ Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALTowerRecPoint * emcRP, AliEMCALDig Int_t iDigit ; - AliEMCALGeometry * geom = gime->EMCALGeometry() ; - 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->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] ; @@ -281,17 +311,33 @@ Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALTowerRecPoint * emcRP, AliEMCALDig //____________________________________________________________________________ void AliEMCALClusterizerv1::GetCalibrationParameters() { - AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; - const AliEMCALDigitizer * dig = gime->Digitizer(BranchName()) ; - - fADCchannelPRE = dig->GetPREchannel() ; - fADCpedestalPRE = dig->GetPREpedestal() ; - - fADCchannelEC = dig->GetECchannel() ; - fADCpedestalEC = dig->GetECpedestal(); - - fADCchannelHC = dig->GetHCchannel() ; - fADCpedestalHC = dig->GetHCpedestal(); + // 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(AliCDBManager::Instance()->IsDefaultStorageSet()){ + AliCDBEntry *entry = (AliCDBEntry*) + AliCDBManager::Instance()->Get("EMCAL/Calib/Data"); + if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject(); + } + 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()); + + fADCchannelECA = dig->GetECAchannel() ; + fADCpedestalECA = dig->GetECApedestal(); + } } //____________________________________________________________________________ @@ -300,504 +346,346 @@ 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 - if ( strcmp(GetTitle(), "") == 0 ) - SetTitle("galice.root") ; + AliRunLoader *rl = AliRunLoader::GetRunLoader(); + fGeom = dynamic_cast(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry(); + fGeom->GetTransformationForSM(); // Global <-> Local + AliInfo(Form("geom 0x%x",fGeom)); - TString branchname = GetName() ; - branchname.Remove(branchname.Index(Version())-1) ; - - AliEMCALGetter * gime = AliEMCALGetter::GetInstance(GetTitle(), branchname.Data(), fToSplit ) ; - if ( gime == 0 ) { - Error("Init", "Could not obtain the Getter object !" ) ; - return ; - } - - fSplitFile = 0 ; - if(fToSplit){ - // construct the name of the file as /path/EMCAL.SDigits.root - //First - extract full path if necessary - TString fileName(GetTitle()) ; - Ssiz_t islash = fileName.Last('/') ; - if(islash(gROOT->GetFile(fileName.Data())); - if(!fSplitFile) - fSplitFile = TFile::Open(fileName.Data(),"update") ; - } - - const AliEMCALGeometry * geom = gime->EMCALGeometry() ; - fNTowers = geom->GetNZ() * geom->GetNPhi() ; if(!gMinuit) gMinuit = new TMinuit(100) ; - gime->PostClusterizer(this) ; - gime->PostRecPoints(branchname ) ; - + fHists = BookHists(); } //____________________________________________________________________________ void AliEMCALClusterizerv1::InitParameters() -{ - fNumberOfPREClusters = fNumberOfECClusters = fNumberOfHCClusters = 0 ; - fPREClusteringThreshold = 0.0001; // must be adjusted according to the noise leve set by digitizer - fECClusteringThreshold = 0.0045; // must be adjusted according to the noise leve set by digitizer - fHCClusteringThreshold = 0.001; // must be adjusted according to the noise leve set by digitizer - fPRELocMaxCut = 0.03 ; - fECLocMaxCut = 0.03 ; - fHCLocMaxCut = 0.03 ; - - fPREW0 = 4.0 ; - fECW0 = 4.5 ; - fHCW0 = 4.5 ; +{ + // Initializes the parameters for the Clusterizer + fNumberOfECAClusters = 0; + fNTowerInGroup = 36; //Produces maximum of 80 pseudoclusters per event + + fECAClusteringThreshold = 0.5; // value obtained from Alexei + fECALocMaxCut = 0.03; // ?? + + fECAW0 = 4.5 ; fTimeGate = 1.e-8 ; - fToUnfold = kFALSE ; - - TString clusterizerName( GetName()) ; - if (clusterizerName.IsNull() ) - clusterizerName = "Default" ; - clusterizerName.Append(":") ; - clusterizerName.Append(Version()) ; - SetName(clusterizerName) ; - fRecPointsInRun = 0 ; + fRecPointsInRun = 0 ; + fMinECut = 0.01; // have to be tune + fCalibData = 0 ; } //____________________________________________________________________________ -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 // = 1 are neighbour - // = 2 are not neighbour but do not continue searching + // = 2 is in different SM, quit from 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::GetInstance()->EMCALGeometry() ; + 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 rowdiff, coldiff; + rv = 0 ; - Int_t rv = 0 ; + fGeom->GetCellIndex(d1->GetId(), nSupMod1,nTower1,nIphi1,nIeta1); + fGeom->GetCellIndex(d2->GetId(), nSupMod2,nTower2,nIphi2,nIeta2); + if(nSupMod1 != nSupMod2) return 2; // different SM - Int_t relid1[4] ; - geom->AbsToRelNumbering(d1->GetId(), relid1) ; + fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nTower1,nIphi1,nIeta1, iphi1,ieta1); + fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nTower2,nIphi2,nIeta2, iphi2,ieta2); - Int_t relid2[4] ; - geom->AbsToRelNumbering(d2->GetId(), relid2) ; + rowdiff = TMath::Abs(iphi1 - iphi2); + coldiff = TMath::Abs(ieta1 - ieta2) ; - if ( (relid1[0] == relid2[0]) && // inside the same EMCAL Arm - (relid1[1]==relid2[1]) ) { // and same tower section - 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() ) < fTimeGate)) - rv = 1 ; - } - else { - if((relid2[2] > relid1[2]) && (relid2[3] > relid1[3]+1)) - rv = 2; // Difference in row numbers is too large to look further - } + if (( coldiff <= 1 ) && ( rowdiff <= 1 )) rv = 1; // neighbours with at least commom vertex + + 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 ; +} - } - else { - - if( (relid1[0] < relid2[0]) || (relid1[1] != relid2[1]) ) - rv=0 ; +//____________________________________________________________________________ +Int_t AliEMCALClusterizerv1::AreInGroup(AliEMCALDigit * d1, AliEMCALDigit * d2) const +{ + // 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, 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 towerGroup1 = -1, towerGroup2 = -1; + rv = 0 ; + + fGeom->GetCellIndex(d1->GetId(), nSupMod1,nTower1,nIphi1,nIeta1); + fGeom->GetCellIndex(d2->GetId(), nSupMod2,nTower2,nIphi2,nIeta2); + if(nSupMod1 != nSupMod2) return 2; // different SM + + static Int_t nTowerInSM = fGeom->GetNCellsInSupMod()/fGeom->GetNCellsInTower(); + + //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(towerGroup1 != towerGroup2) return 3; //different Towergroup - if (gDebug == 2 ) - Info("AreNeighbours", "neighbours=%d, id1=%d, relid1=%d,%d,%d,%d \n id2=%d, relid2=%d,%d,%d,%d", - rv, d1->GetId(), relid1[0], relid1[1], relid1[2], relid1[3], - d2->GetId(), relid2[0], relid2[1], relid2[2], relid2[3]) ; - - return rv ; + //same SM, same towergroup, we're happy + if(towerGroup1 == towerGroup2 && towerGroup2 >= 0) + rv = 1; + + 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); + + 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(Int_t event) +void AliEMCALClusterizerv1::WriteRecPoints() { // Creates new branches with given title // fills and writes into TreeR. - AliEMCALGetter *gime = AliEMCALGetter::GetInstance() ; - - TObjArray * aPRERecPoints = gime->PRERecPoints() ; - TObjArray * aECRecPoints = gime->ECALRecPoints() ; - TObjArray * aHCRecPoints = gime->HCALRecPoints() ; + AliEMCALLoader *emcalLoader = dynamic_cast(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")); - TClonesArray * digits = gime->Digits() ; - TTree * treeR ; - - if(fToSplit){ - if(!fSplitFile) - return ; - fSplitFile->cd() ; - TString name("TreeR") ; - name += event ; - treeR = dynamic_cast(fSplitFile->Get(name)); - } - else{ - treeR = gAlice->TreeR(); - } + TObjArray * aECARecPoints = emcalLoader->RecPoints() ; - if(!treeR){ - gAlice->MakeTree("R", fSplitFile); - treeR = gAlice->TreeR() ; + TClonesArray * digits = emcalLoader->Digits() ; + TTree * treeR = emcalLoader->TreeR(); + if ( treeR==0 ) { + emcalLoader->MakeRecPointsContainer(); + treeR = emcalLoader->TreeR(); } - Int_t index ; - //Evaluate position, dispersion and other RecPoint properties for PRE section - for(index = 0; index < aPRERecPoints->GetEntries(); index++) - (dynamic_cast(aPRERecPoints->At(index)))->EvalAll(fPREW0,digits) ; - aPRERecPoints->Sort() ; - - for(index = 0; index < aPRERecPoints->GetEntries(); index++) - (dynamic_cast(aPRERecPoints->At(index)))->SetIndexInList(index) ; - - aPRERecPoints->Expand(aPRERecPoints->GetEntriesFast()) ; - //Evaluate position, dispersion and other RecPoint properties for EC section - for(index = 0; index < aECRecPoints->GetEntries(); index++) - (dynamic_cast(aECRecPoints->At(index)))->EvalAll(fECW0,digits) ; - - aECRecPoints->Sort() ; - - for(index = 0; index < aECRecPoints->GetEntries(); index++) - (dynamic_cast(aECRecPoints->At(index)))->SetIndexInList(index) ; - - aECRecPoints->Expand(aECRecPoints->GetEntriesFast()) ; - - //Evaluate position, dispersion and other RecPoint properties for HC section - for(index = 0; index < aHCRecPoints->GetEntries(); index++) - (dynamic_cast(aHCRecPoints->At(index)))->EvalAll(fHCW0,digits) ; + for(index = 0; index < aECARecPoints->GetEntries(); index++) + (dynamic_cast(aECARecPoints->At(index)))->EvalAll(fECAW0,digits) ; - aHCRecPoints->Sort() ; + aECARecPoints->Sort() ; - for(index = 0; index < aHCRecPoints->GetEntries(); index++) - (dynamic_cast(aHCRecPoints->At(index)))->SetIndexInList(index) ; + for(index = 0; index < aECARecPoints->GetEntries(); index++) { + (dynamic_cast(aECARecPoints->At(index)))->SetIndexInList(index) ; + (dynamic_cast(aECARecPoints->At(index)))->Print(); + } - aHCRecPoints->Expand(aHCRecPoints->GetEntriesFast()) ; - Int_t bufferSize = 32000 ; Int_t splitlevel = 0 ; - - //PRE section branch - TBranch * branchPRE = treeR->Branch("EMCALPRERP","TObjArray",&aPRERecPoints,bufferSize,splitlevel); - branchPRE->SetTitle(BranchName()); //EC section branch - TBranch * branchEC = treeR->Branch("EMCALECRP","TObjArray",&aECRecPoints,bufferSize,splitlevel); - branchEC->SetTitle(BranchName()); + + TBranch * branchECA = 0; + if ((branchECA = treeR->GetBranch("EMCALECARP"))) + branchECA->SetAddress(&aECARecPoints); + else + treeR->Branch("EMCALECARP","TObjArray",&aECARecPoints,bufferSize,splitlevel); + treeR->Fill() ; + + emcalLoader->WriteRecPoints("OVERWRITE"); - //HC section branch - TBranch * branchHC = treeR->Branch("EMCALHCRP","TObjArray",&aHCRecPoints,bufferSize,splitlevel); - branchHC->SetTitle(BranchName()); - - //And Finally clusterizer branch - AliEMCALClusterizerv1 * cl = (AliEMCALClusterizerv1*)gime->Clusterizer(BranchName()) ; - TBranch * clusterizerBranch = treeR->Branch("AliEMCALClusterizer","AliEMCALClusterizerv1", - &cl,bufferSize,splitlevel); - clusterizerBranch->SetTitle(BranchName()); - - branchPRE->Fill() ; - branchEC->Fill() ; - branchHC->Fill() ; - clusterizerBranch->Fill() ; - - treeR->AutoSave() ; //Write(0,kOverwrite) ; - if(gAlice->TreeR()!=treeR) - treeR->Delete(); } //____________________________________________________________________________ -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 + + AliEMCALLoader *emcalLoader = dynamic_cast(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")); + + TObjArray * aECARecPoints = emcalLoader->RecPoints() ; - AliEMCALGetter * gime = AliEMCALGetter::GetInstance() ; + if (fGeom==0) + AliFatal("Did not get geometry from EMCALLoader"); - AliEMCALGeometry * geom = gime->EMCALGeometry() ; + aECARecPoints->Clear(); + //Start with pseudoclusters, if option + if(strstr(option,"pseudo")) { + TClonesArray * digits = emcalLoader->Digits() ; + TClonesArray * digitsC = dynamic_cast(digits->Clone()); - TObjArray * aPRERecPoints = gime->PRERecPoints(BranchName()) ; - TObjArray * aECRecPoints = gime->ECALRecPoints(BranchName()) ; - TObjArray * aHCRecPoints = gime->HCALRecPoints(BranchName()) ; + TIter nextdigit(digitsC) ; + AliEMCALDigit * digit; - aPRERecPoints->Delete() ; - aECRecPoints->Delete() ; - aHCRecPoints->Delete() ; + AliEMCALRecPoint * recPoint = 0 ; + int ineb=0; + nextdigit.Reset(); - TClonesArray * digits = gime->Digits() ; - if ( !digits ) { - Fatal("MakeClusters -> Digits with name %s not found", GetName() ) ; - } + // PseudoClusterization starts + while ( (digit = dynamic_cast(nextdigit())) ) { // scan over the list of digitsC + recPoint = 0 ; + TArrayI clusterECAdigitslist(1000); // what is this - TIter next(digits) ; - AliEMCALDigit * digit ; - Int_t ndigEC=0, ndigPRE=0, ndigHC=0 ; - - // count the number of digits in EC section - while ( (digit = dynamic_cast(next())) ) { // scan over the list of digits - if (geom->IsInECAL(digit->GetId())) - ndigEC++ ; - else if (geom->IsInPRE(digit->GetId())) - ndigPRE++; - else if (geom->IsInHCAL(digit->GetId())) - ndigHC++; - else { - Error("MakeClusters", "id = %d is a wrong ID!", digit->GetId()) ; - abort() ; - } - } + 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++ ; - // add amplitude of PRE and EC sections - Int_t digEC ; - for (digEC = 0 ; digEC < ndigEC ; digEC++) { - digit = dynamic_cast(digits->At(digEC)) ; - Int_t digPRE ; - for (digPRE = ndigEC ; digPRE < ndigEC+ndigPRE ; digPRE++) { - AliEMCALDigit * digitPRE = dynamic_cast(digits->At(digPRE)) ; - if ( geom->AreInSameTower(digit->GetId(), digitPRE->GetId()) ){ - Float_t amp = static_cast(digit->GetAmp()) + geom->GetSummationFraction() * static_cast(digitPRE->GetAmp()) + 0.5 ; - digit->SetAmp(static_cast(amp)) ; - break ; + recPoint->SetClusterType(AliESDCaloCluster::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 } } - if (gDebug) - Info("MakeClusters","id = %d amp = %d", digit->GetId(), digit->GetAmp()) ; - } + if(recPoint) + AliDebug(1,Form("MakeClusters: cl.e %f \n", recPoint->GetEnergy())); + //if(recPoint) cout << "cl.e " << recPoint->GetEnergy() << endl; + delete digitsC ; + } - TClonesArray * digitsC = dynamic_cast(digits->Clone()) ; - - - // Clusterization starts - - TIter nextdigit(digitsC) ; - Bool_t notremovedEC = kTRUE, notremovedPRE = kTRUE ; - - while ( (digit = dynamic_cast(nextdigit())) ) { // scan over the list of digitsC - AliEMCALRecPoint * clu = 0 ; - - TArrayI clusterPREdigitslist(50), clusterECdigitslist(50), clusterHCdigitslist(50); - - Bool_t inPRE = kFALSE, inECAL = kFALSE, inHCAL = kFALSE ; - if( geom->IsInPRE(digit->GetId()) ) { - inPRE = kTRUE ; - } - else if( geom->IsInECAL(digit->GetId()) ) { - inECAL = kTRUE ; - } - else if( geom->IsInHCAL(digit->GetId()) ) { - inHCAL = kTRUE ; - } - - if (gDebug == 2) { - if (inPRE) - Info("MakeClusters","id = %d, ene = %f , thre = %f ", - digit->GetId(),Calibrate(digit->GetAmp(), 0), fPREClusteringThreshold) ; - if (inECAL) - Info("MakeClusters","id = %d, ene = %f , thre = %f", - digit->GetId(),Calibrate(digit->GetAmp(), 1), fECClusteringThreshold) ; - if (inHCAL) - Info("MakeClusters","id = %d, ene = %f , thre = %f", - digit->GetId(),Calibrate(digit->GetAmp(), 2), fHCClusteringThreshold ) ; - } - - if ( (inPRE && (Calibrate(digit->GetAmp(), 0) > fPREClusteringThreshold )) || - (inECAL && (Calibrate(digit->GetAmp(), 1) > fECClusteringThreshold )) || - (inHCAL && (Calibrate(digit->GetAmp(), 2) > fHCClusteringThreshold )) ) { - - Int_t iDigitInPRECluster = 0, iDigitInECCluster = 0, iDigitInHCCluster = 0; - Int_t where ; // PRE = 0, ECAl = 1, HCAL = 2 - - // Find the seed in each of the section ECAL/PRE/HCAL - - if( geom->IsInECAL(digit->GetId()) ) { - where = 1 ; // to tell we are in ECAL - // start a new Tower RecPoint - if(fNumberOfECClusters >= aECRecPoints->GetSize()) - aECRecPoints->Expand(2*fNumberOfECClusters+1) ; - AliEMCALTowerRecPoint * rp = new AliEMCALTowerRecPoint("") ; - rp->SetECAL() ; - aECRecPoints->AddAt(rp, fNumberOfECClusters) ; - clu = dynamic_cast(aECRecPoints->At(fNumberOfECClusters)) ; - fNumberOfECClusters++ ; - clu->AddDigit(*digit, Calibrate(digit->GetAmp(), where)) ; - clusterECdigitslist[iDigitInECCluster] = digit->GetIndexInList() ; - iDigitInECCluster++ ; - digitsC->Remove(digit) ; - if (gDebug == 2 ) - Info("MakeClusters","OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp(), 1), fECClusteringThreshold) ; - - } - else if( geom->IsInPRE(digit->GetId()) ) { - where = 0 ; // to tell we are in PRE - // start a new Pre Shower cluster - if(fNumberOfPREClusters >= aPRERecPoints->GetSize()) - aPRERecPoints->Expand(2*fNumberOfPREClusters+1); - AliEMCALTowerRecPoint * rp = new AliEMCALTowerRecPoint("") ; - rp->SetPRE() ; - aPRERecPoints->AddAt(rp, fNumberOfPREClusters) ; - clu = dynamic_cast(aPRERecPoints->At(fNumberOfPREClusters)) ; - fNumberOfPREClusters++ ; - clu->AddDigit(*digit, Calibrate(digit->GetAmp(), where)); - clusterPREdigitslist[iDigitInPRECluster] = digit->GetIndexInList() ; - iDigitInPRECluster++ ; - digitsC->Remove(digit) ; - if (gDebug == 2 ) - Info("MakeClusters","OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp(), 0), fPREClusteringThreshold) ; - - nextdigit.Reset() ; - - // Here we remove remaining EC digits, which cannot make a cluster - - if( notremovedEC ) { - while( ( digit = dynamic_cast(nextdigit()) ) ) { - if( geom->IsInECAL(digit->GetId()) ) - digitsC->Remove(digit) ; - else - break ; - } - notremovedEC = kFALSE ; - } + //Now do real clusters + TClonesArray * digits = emcalLoader->Digits() ; + TClonesArray * digitsC = dynamic_cast(digits->Clone()); - } - else if( geom->IsInHCAL(digit->GetId()) ) { - where = 2 ; // to tell we are in HCAL - // start a new HCAL cluster - if(fNumberOfHCClusters >= aHCRecPoints->GetSize()) - aHCRecPoints->Expand(2*fNumberOfHCClusters+1); - AliEMCALTowerRecPoint * rp = new AliEMCALTowerRecPoint("") ; - rp->SetHCAL() ; - aHCRecPoints->AddAt(rp, fNumberOfHCClusters) ; - clu = dynamic_cast(aHCRecPoints->At(fNumberOfHCClusters)) ; - fNumberOfHCClusters++ ; - clu->AddDigit(*digit, Calibrate(digit->GetAmp(), where)); - clusterHCdigitslist[iDigitInHCCluster] = digit->GetIndexInList() ; - iDigitInHCCluster++ ; - digitsC->Remove(digit) ; - if (gDebug == 2 ) - Info("MakeClusters","OK id = %d, ene = %f , thre = %f ", digit->GetId(),Calibrate(digit->GetAmp(), 2), fHCClusteringThreshold) ; - - nextdigit.Reset() ; - - // Here we remove remaining PRE digits, which cannot make a cluster - - if( notremovedPRE ) { - while( ( digit = dynamic_cast(nextdigit()) ) ) { - if( geom->IsInPRE(digit->GetId()) ) - digitsC->Remove(digit) ; - else - break ; - } - notremovedPRE = kFALSE ; - } - } - - nextdigit.Reset() ; - - AliEMCALDigit * digitN ; - Int_t index = 0 ; + TIter nextdigit(digitsC) ; + AliEMCALDigit * digit; - // Do the Clustering in each of the three section ECAL/PRE/HCAL + AliEMCALRecPoint * recPoint = 0 ; + int ineb=0; + nextdigit.Reset(); - while (index < iDigitInECCluster){ // scan over digits already in cluster - digit = (AliEMCALDigit*)digits->At(clusterECdigitslist[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 !!!!! - // Info("MakeClusters","id1 = %d, id2 = %d , neighbours = %d", digit->GetId(), digitN->GetId(), ineb) ; - switch (ineb ) { - case 0 : // not a neighbour - break ; - case 1 : // are neighbours - clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), 1) ) ; - clusterECdigitslist[iDigitInECCluster] = digitN->GetIndexInList() ; - iDigitInECCluster++ ; - digitsC->Remove(digitN) ; - break ; - case 2 : // too far from each other - goto endofloop1; - } // switch - - } // while digitN - - endofloop1: ; - nextdigit.Reset() ; - } // loop over EC cluster + double e=0.0, ehs = 0.0; + while ( (digit = dynamic_cast(nextdigit())) ) { // 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; + } + AliDebug(1,Form("MakeClusters: Number of digits %d -> (e %f), ehs %d\n", + digits->GetEntries(),fMinECut,ehs)); + //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(fGeom->CheckAbsCellId(digit->GetId()) && (Calibrate(digit->GetAmp(), digit->GetId()) > fECAClusteringThreshold ) ){ + 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(AliESDCaloCluster::kClusterv1); + + recPoint->AddDigit(*digit, Calibrate(digit->GetAmp(), digit->GetId())) ; + clusterECAdigitslist[iDigitInECACluster] = digit->GetIndexInList() ; + iDigitInECACluster++ ; + digitsC->Remove(digit) ; + AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , thre = %f \n", digit->GetId(), + Calibrate(digit->GetAmp(),digit->GetId()), fECAClusteringThreshold)); + nextdigit.Reset(); // will start from beggining - index = 0 ; - while (index < iDigitInPRECluster){ // scan over digits already in cluster - digit = (AliEMCALDigit*)digits->At(clusterPREdigitslist[index]) ; + 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 - Int_t ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!! - // Info("MakeClusters","id1 = %d, id2 = %d , neighbours = %d", digit->GetId(), digitN->GetId(), ineb) ; - switch (ineb ) { - case 0 : // not a neighbour - break ; - case 1 : // are neighbours - clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), 0) ) ; - clusterPREdigitslist[iDigitInPRECluster] = digitN->GetIndexInList() ; - iDigitInPRECluster++ ; - digitsC->Remove(digitN) ; - break ; - case 2 : // too far from each other - goto endofloop2; + 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 - - } // while digitN - - endofloop2: ; - nextdigit.Reset() ; - } // loop over PRE cluster - - index = 0 ; - while (index < iDigitInHCCluster){ // scan over digits already in cluster - digit = (AliEMCALDigit*)digits->At(clusterHCdigitslist[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 !!!!! - //Info("MakeClusters","id1 = %d, id2 = %d , neighbours = %d", digit->GetId(), digitN->GetId(), ineb) ; - switch (ineb ) { - case 0 : // not a neighbour - break ; - case 1 : // are neighbours - clu->AddDigit(*digitN, Calibrate( digitN->GetAmp(), 2) ) ; - clusterHCdigitslist[iDigitInHCCluster] = digitN->GetIndexInList() ; - iDigitInHCCluster++ ; - digitsC->Remove(digitN) ; - break ; - case 2 : // too far from each other - goto endofloop3; - } // switch - } // while digitN - - endofloop3: ; - nextdigit.Reset() ; - } // loop over HC cluster - - } // energy theshold - } // while digit + } // scan over the reduced list of digits + } // scan over digits already in cluster + nextdigit.Reset() ; // will start from beggining + } + } // while digit + if(recPoint) + AliDebug(1,Form("MakeClusters: cl.e %f \n", recPoint->GetEnergy())); + //if(recPoint) cout << "cl.e " << recPoint->GetEnergy() << endl; delete digitsC ; + + AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,digits->GetEntriesFast())); } //____________________________________________________________________________ -void AliEMCALClusterizerv1::MakeUnfolding() +void AliEMCALClusterizerv1::MakeUnfolding() const { Fatal("AliEMCALClusterizerv1::MakeUnfolding", "--> Unfolding not implemented") ; @@ -816,10 +704,10 @@ Double_t AliEMCALClusterizerv1::ShowerShape(Double_t r) } //____________________________________________________________________________ -void AliEMCALClusterizerv1::UnfoldCluster(AliEMCALTowerRecPoint * iniTower, - Int_t nMax, - AliEMCALDigit ** maxAt, - Float_t * maxAtEnergy) +void AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * /*iniTower*/, + Int_t /*nMax*/, + AliEMCALDigit ** /*maxAt*/, + Float_t * /*maxAtEnergy*/) const { // Performs the unfolding of a cluster with nMax overlapping showers @@ -828,7 +716,9 @@ void AliEMCALClusterizerv1::UnfoldCluster(AliEMCALTowerRecPoint * iniTower, } //_____________________________________________________________________________ -void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Double_t & fret, Double_t * x, Int_t iflag) +void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & /*nPar*/, Double_t * /*Grad*/, + Double_t & /*fret*/, + Double_t * /*x*/, Int_t /*iflag*/) { // Calculates the Chi square for the cluster unfolding minimization // Number of parameters, Gradient, Chi squared, parameters, what to do @@ -836,7 +726,7 @@ void AliEMCALClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Do ::Fatal("UnfoldingChiSquare","Unfolding not implemented") ; } //____________________________________________________________________________ -void AliEMCALClusterizerv1::Print(Option_t * option)const +void AliEMCALClusterizerv1::Print(Option_t * /*option*/)const { // Print clusterizer parameters @@ -846,106 +736,60 @@ void AliEMCALClusterizerv1::Print(Option_t * option)const // Print parameters - TString taskName(GetName()) ; + TString taskName(GetName()) ; taskName.ReplaceAll(Version(), "") ; - message += "--------------- " ; - message += taskName.Data() ; - message += " " ; - message += GetTitle() ; - message += "-----------\n" ; - message += "Clusterizing digits from the file: " ; - message += taskName.Data() ; - message += "\n Branch: " ; - message += GetName() ; - message += "\n Pre Shower Clustering threshold = " ; - message += fPREClusteringThreshold ; - message += "\n Pre Shower Local Maximum cut = " ; - message += fPRELocMaxCut ; - message += "\n Pre Shower Logarothmic weight = " ; - message += fPREW0 ; - message += "\n EC Clustering threshold = " ; - message += fECClusteringThreshold ; - message += "\n EC Local Maximum cut = " ; - message += fECLocMaxCut ; - message += "\n EC Logarothmic weight = " ; - message += fECW0 ; - message += "\n Pre Shower Clustering threshold = " ; - message += fHCClusteringThreshold ; - message += "\n HC Local Maximum cut = " ; - message += fHCLocMaxCut ; - message += "\n HC Logarothmic weight = " ; - message += fHCW0 ; + printf("--------------- "); + printf(taskName.Data()) ; + printf(" "); + printf(GetTitle()) ; + printf("-----------\n"); + printf("Clusterizing digits from the file: "); + printf(taskName.Data()); + printf("\n Branch: "); + printf(GetName()); + printf("\n ECA Local Maximum cut = %f", fECALocMaxCut); + printf("\n ECA Logarithmic weight = %f", fECAW0); if(fToUnfold) - message +="\nUnfolding on\n" ; + printf("\nUnfolding on\n"); else - message += "\nUnfolding off\n"; + printf("\nUnfolding off\n"); - message += "------------------------------------------------------------------" ; + printf("------------------------------------------------------------------"); } else - message += "AliEMCALClusterizerv1 not initialized " ; - - Info("Print", message.Data() ) ; + printf("AliEMCALClusterizerv1 not initialized ") ; } //____________________________________________________________________________ void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option) { // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer - - TObjArray * aPRERecPoints = AliEMCALGetter::GetInstance()->PRERecPoints() ; - TObjArray * aECRecPoints = AliEMCALGetter::GetInstance()->ECALRecPoints() ; - TObjArray * aHCRecPoints = AliEMCALGetter::GetInstance()->HCALRecPoints() ; - - Info("PrintRecPoints", "Clusterization result:") ; + AliEMCALLoader *emcalLoader = dynamic_cast(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")); + TObjArray * aECARecPoints = emcalLoader->RecPoints() ; + + printf("PrintRecPoints: Clusterization result:") ; - printf("event # %d\n", gAlice->GetEvNumber() ) ; - printf(" Found %d PRE SHOWER RecPoints, %d EC Rec Points and %d HC Rec Points\n ", - aPRERecPoints->GetEntriesFast(), aECRecPoints->GetEntriesFast(), aHCRecPoints->GetEntriesFast() ) ; + printf("event # %d\n", emcalLoader->GetRunLoader()->GetEventNumber() ) ; + printf(" Found %d ECA Rec Points\n ", + aECARecPoints->GetEntriesFast()) ; - fRecPointsInRun += aPRERecPoints->GetEntriesFast() ; - fRecPointsInRun += aECRecPoints->GetEntriesFast() ; - fRecPointsInRun += aHCRecPoints->GetEntriesFast() ; + fRecPointsInRun += aECARecPoints->GetEntriesFast() ; if(strstr(option,"all")) { - - //Pre shower recPoints - - printf("-----------------------------------------------------------------------\n") ; - printf("Clusters in PRE section\n") ; - printf("Index Ene(GeV) Multi Module phi r theta X Y Z Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ; - - Int_t index ; - - for (index = 0 ; index < aPRERecPoints->GetEntries() ; index++) { - AliEMCALTowerRecPoint * rp = dynamic_cast(aPRERecPoints->At(index)) ; - TVector3 globalpos; - rp->GetGlobalPosition(globalpos); - TVector3 localpos; - rp->GetLocalPosition(localpos); - Float_t lambda[2]; - rp->GetElipsAxis(lambda); - Int_t * primaries; - Int_t nprimaries; - primaries = rp->GetPrimaries(nprimaries); - printf("\n%6d %8.4f %3d %2d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ", - rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetEMCALArm(), - globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(), - rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ; - for (Int_t iprimary=0; iprimaryGetEntries() ; index++) { - AliEMCALTowerRecPoint * rp = dynamic_cast(aECRecPoints->At(index)) ; + Float_t maxE=0; + Float_t maxL1=0; + Float_t maxL2=0; + Float_t maxDis=0; + + for (index = 0 ; index < aECARecPoints->GetEntries() ; index++) { + AliEMCALRecPoint * rp = dynamic_cast(aECARecPoints->At(index)) ; TVector3 globalpos; - rp->GetGlobalPosition(globalpos); + //rp->GetGlobalPosition(globalpos); TVector3 localpos; rp->GetLocalPosition(localpos); Float_t lambda[2]; @@ -953,39 +797,67 @@ void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option) Int_t * primaries; Int_t nprimaries; primaries = rp->GetPrimaries(nprimaries); - printf("\n%6d %8.4f %3d %2d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ", - rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetEMCALArm(), + 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) ; + ///////////// + if(rp->GetEnergy()>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()); + ///////////// for (Int_t iprimary=0; iprimaryGetEntries() ; index++) { - AliEMCALTowerRecPoint * rp = dynamic_cast(aHCRecPoints->At(index)) ; - TVector3 globalpos; - rp->GetGlobalPosition(globalpos); - TVector3 localpos; - rp->GetLocalPosition(localpos); - Float_t lambda[2]; - rp->GetElipsAxis(lambda); - Int_t * primaries; - Int_t nprimaries; - primaries = rp->GetPrimaries(nprimaries); - printf("\n%6d %8.4f %3d %2d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ", - rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), rp->GetEMCALArm(), - globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(), - rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ; - for (Int_t iprimary=0; iprimaryFill(maxE); + fMaxL1->Fill(maxL1); + fMaxL2->Fill(maxL2); + fMaxDis->Fill(maxDis); + 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 + // + 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 + + return AliEMCALHistoUtilities::MoveHistsToList("EmcalClusterizerv1ControlHists", kFALSE); +} + +void AliEMCALClusterizerv1::SaveHists(const char *fn) +{ + AliEMCALHistoUtilities::SaveListOfHists(fHists, fn, kTRUE); +} + +void AliEMCALClusterizerv1::Browse(TBrowser* b) +{ + if(fHists) b->Add(fHists); + TTask::Browse(b); +}