X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliEMCALClusterizerv1.cxx;h=700aacc0857cdd5303d3247f992dad84bc3201de;hb=8221b41b8b59a093459af5fb8dc493ce12117bbd;hp=66186ebcc8210c85dbf4f18179746fbf36321686;hpb=fb43ada4db956ceb59baa593557bff7082e4bfc2;p=u%2Fmrichter%2FAliRoot.git diff --git a/EMCAL/AliEMCALClusterizerv1.cxx b/EMCAL/AliEMCALClusterizerv1.cxx index 66186ebcc82..700aacc0857 100644 --- a/EMCAL/AliEMCALClusterizerv1.cxx +++ b/EMCAL/AliEMCALClusterizerv1.cxx @@ -46,69 +46,165 @@ // --- 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 "AliEMCALGetter.h" +#include "AliRunLoader.h" +#include "AliRun.h" +#include "AliESD.h" +#include "AliEMCALLoader.h" #include "AliEMCALClusterizerv1.h" -#include "AliEMCALTowerRecPoint.h" +#include "AliEMCALRecPoint.h" #include "AliEMCALDigit.h" #include "AliEMCALDigitizer.h" #include "AliEMCAL.h" #include "AliEMCALGeometry.h" +#include "AliEMCALHistoUtilities.h" +#include "AliCDBManager.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),fTimeGate(0.),fMinECut(0.) { // default ctor (to be used mainly by Streamer) InitParameters() ; - fDefaultInit = kTRUE ; + fGeom = AliEMCALGeometry::GetInstance(); + 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),fTimeGate(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), + fTimeGate(clus.fTimeGate), + fMinECut(clus.fMinECut) +{ + //copy ctor } //____________________________________________________________________________ AliEMCALClusterizerv1::~AliEMCALClusterizerv1() { // dtor - } //____________________________________________________________________________ const TString AliEMCALClusterizerv1::BranchName() const { return GetName(); - } //____________________________________________________________________________ -Float_t AliEMCALClusterizerv1::Calibrate(Int_t amp) const +Float_t AliEMCALClusterizerv1::Calibrate(Int_t amp, Int_t AbsId) { - //To be replased later by the method, reading individual parameters from the database - return -fADCpedestalECA + amp * fADCchannelECA ; + + // 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 ; + } //____________________________________________________________________________ @@ -125,23 +221,27 @@ void AliEMCALClusterizerv1::Exec(Option_t * option) if(strstr(option,"print")) Print("") ; - AliEMCALGetter * gime = AliEMCALGetter::Instance(GetTitle()) ; + 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 = gime->MaxEvent() - 1 ; - else - fLastEvent = TMath::Min(fFirstEvent, gime->MaxEvent()); + fLastEvent = rl->GetNumberOfEvents() - 1; Int_t nEvents = fLastEvent - fFirstEvent + 1; Int_t ievent ; + rl->LoadDigits("EMCAL"); for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) { - gime->Event(ievent,"D") ; - - GetCalibrationParameters() ; + rl->GetEvent(ievent); fNumberOfECAClusters = 0; - - MakeClusters() ; + + if(strstr(option,"pseudo")) + MakeClusters("pseudo") ; //both types + else + MakeClusters("") ; //only the real clusters if(fToUnfold) MakeUnfolding() ; @@ -152,7 +252,7 @@ void AliEMCALClusterizerv1::Exec(Option_t * option) PrintRecPoints(option) ; //increment the total number of recpoints per run - fRecPointsInRun += gime->ECARecPoints()->GetEntriesFast() ; + fRecPointsInRun += emcalLoader->RecPoints()->GetEntriesFast() ; } Unload(); @@ -160,21 +260,22 @@ void AliEMCALClusterizerv1::Exec(Option_t * option) 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 ) ; + 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::Instance() ; - 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 gMinuit->SetFCN(AliEMCALClusterizerv1::UnfoldingChiSquare) ; @@ -194,16 +295,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[3] ; 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] ; @@ -258,15 +358,33 @@ Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALTowerRecPoint * emcRP, AliEMCALDig //____________________________________________________________________________ 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(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(); + } } //____________________________________________________________________________ @@ -275,16 +393,15 @@ 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()); - - AliEMCALGeometry * geom = gime->EMCALGeometry() ; + AliRunLoader *rl = AliRunLoader::GetRunLoader(); + fGeom = dynamic_cast(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry(); + fGeom->GetTransformationForSM(); // Global <-> Local + AliInfo(Form("geom 0x%x",fGeom)); - fNTowers = geom->GetNZ() * geom->GetNPhi() ; if(!gMinuit) gMinuit = new TMinuit(100) ; - if ( !gime->Clusterizer() ) - gime->PostClusterizer(this); + fHists = BookHists(); } //____________________________________________________________________________ @@ -292,69 +409,112 @@ 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 ; + + 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 ; 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::Instance()->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[3] ; - geom->AbsToRelNumbering(d1->GetId(), relid1) ; + fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nTower1,nIphi1,nIeta1, iphi1,ieta1); + fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nTower2,nIphi2,nIeta2, iphi2,ieta2); - Int_t relid2[3] ; - 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 - Int_t rowdiff = TMath::Abs( relid1[1] - relid2[1] ) ; - Int_t coldiff = TMath::Abs( relid1[2] - relid2[2] ) ; - - if (( coldiff <= 1 ) && ( rowdiff <= 1 )){ - if(TMath::Abs(d1->GetTime() - d2->GetTime() ) < fTimeGate) - rv = 1 ; - } - else { - if((relid2[1] > relid1[1]) && (relid2[2] > relid1[2]+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]) - 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 ) - printf("AreNeighbours: neighbours=%d, id1=%d, relid1=%d,%d,%d \n id2=%d, relid2=%d,%d,%d ", - rv, d1->GetId(), relid1[0], relid1[1], relid1[2], - d2->GetId(), relid2[0], relid2[1], relid2[2]) ; - - 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 - AliEMCALGetter * gime = AliEMCALGetter::Instance() ; - gime->EmcalLoader()->UnloadDigits() ; - gime->EmcalLoader()->UnloadRecPoints() ; + AliEMCALLoader *emcalLoader = dynamic_cast(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")); + + emcalLoader->UnloadDigits() ; + emcalLoader->UnloadRecPoints() ; } //____________________________________________________________________________ @@ -364,141 +524,211 @@ void AliEMCALClusterizerv1::WriteRecPoints() // Creates new branches with given title // fills and writes into TreeR. - AliEMCALGetter *gime = AliEMCALGetter::Instance() ; + AliEMCALLoader *emcalLoader = dynamic_cast(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")); - TObjArray * aECARecPoints = gime->ECARecPoints() ; + TObjArray * aECARecPoints = emcalLoader->RecPoints() ; - TClonesArray * digits = gime->Digits() ; - TTree * treeR = gime->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 EC section for(index = 0; index < aECARecPoints->GetEntries(); index++) - (dynamic_cast(aECARecPoints->At(index)))->EvalAll(fECAW0,digits) ; + (dynamic_cast(aECARecPoints->At(index)))->EvalAll(fECAW0,digits) ; aECARecPoints->Sort() ; - for(index = 0; index < aECARecPoints->GetEntries(); index++) - (dynamic_cast(aECARecPoints->At(index)))->SetIndexInList(index) ; + for(index = 0; index < aECARecPoints->GetEntries(); index++) { + (dynamic_cast(aECARecPoints->At(index)))->SetIndexInList(index) ; + (dynamic_cast(aECARecPoints->At(index)))->Print(); + } - aECARecPoints->Expand(aECARecPoints->GetEntriesFast()) ; - Int_t bufferSize = 32000 ; Int_t splitlevel = 0 ; //EC section branch - TBranch * branchECA = treeR->Branch("EMCALECARP","TObjArray",&aECARecPoints,bufferSize,splitlevel); - branchECA->SetTitle(BranchName()); + + TBranch * branchECA = 0; + if ((branchECA = treeR->GetBranch("EMCALECARP"))) + branchECA->SetAddress(&aECARecPoints); + else + treeR->Branch("EMCALECARP","TObjArray",&aECARecPoints,bufferSize,splitlevel); + treeR->Fill() ; - branchECA->Fill() ; + emcalLoader->WriteRecPoints("OVERWRITE"); - gime->WriteRecPoints("OVERWRITE"); - gime->WriteClusterizer("OVERWRITE"); } //____________________________________________________________________________ -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::Instance() ; + if (fGeom==0) + AliFatal("Did not get geometry from EMCALLoader"); - AliEMCALGeometry * geom = gime->EMCALGeometry() ; + aECARecPoints->Clear(); - TObjArray * aECARecPoints = gime->ECARecPoints() ; + //Start with pseudoclusters, if option + if(strstr(option,"pseudo")) { + TClonesArray * digits = emcalLoader->Digits() ; + TClonesArray * digitsC = dynamic_cast(digits->Clone()); - aECARecPoints->Delete() ; + TIter nextdigit(digitsC) ; + AliEMCALDigit * digit; - TClonesArray * digits = gime->Digits() ; + AliEMCALRecPoint * recPoint = 0 ; + int ineb=0; + nextdigit.Reset(); - TIter next(digits) ; - AliEMCALDigit * digit ; - Int_t ndigECA=0 ; + // PseudoClusterization starts + while ( (digit = dynamic_cast(nextdigit())) ) { // scan over the list of digitsC + recPoint = 0 ; + TArrayI clusterECAdigitslist(1000); // what is this - // 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() ; - } - } - 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) ; - } - 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) ; - AliEMCALTowerRecPoint * rp = new AliEMCALTowerRecPoint("") ; - rp->SetECA() ; + 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) ; - clu = dynamic_cast(aECARecPoints->At(fNumberOfECAClusters)) ; - fNumberOfECAClusters++ ; - clu->AddDigit(*digit, Calibrate(digit->GetAmp())) ; + recPoint = dynamic_cast(aECARecPoints->At(fNumberOfECAClusters)) ; + fNumberOfECAClusters++ ; + + recPoint->SetClusterType(AliESDCaloCluster::kPseudoCluster); + + recPoint->AddDigit(*digit, 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("MakePseudoClusters: OK id = %d, adc = %f \n", digit->GetId(), digit->GetAmp())); + nextdigit.Reset(); // will start from beggining - AliEMCALDigit * digitN ; - Int_t index = 0 ; + 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(recPoint) + AliDebug(1,Form("MakeClusters: cl.e %f \n", recPoint->GetEnergy())); + //if(recPoint) cout << "cl.e " << recPoint->GetEnergy() << endl; + delete digitsC ; + } - // Do the Clustering + //Now do real clusters + TClonesArray * digits = emcalLoader->Digits() ; + TClonesArray * digitsC = dynamic_cast(digits->Clone()); + TIter nextdigit(digitsC) ; + AliEMCALDigit * digit; + + AliEMCALRecPoint * recPoint = 0 ; + int ineb=0; + nextdigit.Reset(); + + 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 + + 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]) ; + 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; + 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 - - endofloop1: ; - nextdigit.Reset() ; - } // loop over ECA 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())); } //____________________________________________________________________________ @@ -521,7 +751,7 @@ Double_t AliEMCALClusterizerv1::ShowerShape(Double_t r) } //____________________________________________________________________________ -void AliEMCALClusterizerv1::UnfoldCluster(AliEMCALTowerRecPoint * /*iniTower*/, +void AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * /*iniTower*/, Int_t /*nMax*/, AliEMCALDigit ** /*maxAt*/, Float_t * /*maxAtEnergy*/) const @@ -582,11 +812,12 @@ 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() ; + AliEMCALLoader *emcalLoader = dynamic_cast(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")); + TObjArray * aECARecPoints = emcalLoader->RecPoints() ; + printf("PrintRecPoints: Clusterization result:") ; - printf("event # %d\n", gAlice->GetEvNumber() ) ; + printf("event # %d\n", emcalLoader->GetRunLoader()->GetEventNumber() ) ; printf(" Found %d ECA Rec Points\n ", aECARecPoints->GetEntriesFast()) ; @@ -597,11 +828,15 @@ void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option) 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") ; - + Float_t maxE=0; + Float_t maxL1=0; + Float_t maxL2=0; + Float_t maxDis=0; + for (index = 0 ; index < aECARecPoints->GetEntries() ; index++) { - AliEMCALTowerRecPoint * rp = dynamic_cast(aECARecPoints->At(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]; @@ -609,14 +844,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; 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); +}