X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliEMCALClusterizerv1.cxx;h=9cd9bcdda9309a3a8a8b3c626e561ed76faf9660;hb=dd70110949cdcdc8c7c1cba590a9a2ea7e6de638;hp=0085172be7352765c07c0260a7f03190c94ec446;hpb=d87bd0451ba0a37ae17dc2015f4d8dcb91a6db90;p=u%2Fmrichter%2FAliRoot.git diff --git a/EMCAL/AliEMCALClusterizerv1.cxx b/EMCAL/AliEMCALClusterizerv1.cxx index 0085172be73..9cd9bcdda93 100644 --- a/EMCAL/AliEMCALClusterizerv1.cxx +++ b/EMCAL/AliEMCALClusterizerv1.cxx @@ -1,4 +1,3 @@ - /************************************************************************** * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * * * @@ -19,6 +18,7 @@ //*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute) // August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction // of new IO (à la PHOS) +// Mar 2007, Aleksei Pavlinov - new algoritmh of pseudo clusters ////////////////////////////////////////////////////////////////////////////// // Clusterization class. Performs clusterization (collects neighbouring active cells) and // unfolds the clusters having several local maxima. @@ -46,145 +46,208 @@ // // 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 "AliEMCALLoader.h" #include "AliEMCALClusterizerv1.h" #include "AliEMCALRecPoint.h" #include "AliEMCALDigit.h" #include "AliEMCALDigitizer.h" #include "AliEMCAL.h" #include "AliEMCALGeometry.h" +#include "AliEMCALRawUtils.h" +#include "AliEMCALHistoUtilities.h" +#include "AliCDBManager.h" -ClassImp(AliEMCALClusterizerv1) +class AliCDBStorage; +#include "AliCDBEntry.h" -Int_t addOn[20][60][60]; +ClassImp(AliEMCALClusterizerv1) //____________________________________________________________________________ - AliEMCALClusterizerv1::AliEMCALClusterizerv1() : AliEMCALClusterizer() +AliEMCALClusterizerv1::AliEMCALClusterizerv1() + : AliEMCALClusterizer(), + fHists(0),fPointE(0),fPointL1(0),fPointL2(0), + fPointDis(0),fPointMult(0),fDigitAmp(0),fMaxE(0), + fMaxL1(0),fMaxL2(0),fMaxDis(0),fGeom(0), + fDefaultInit(kTRUE), + fToUnfold(kFALSE), + fNumberOfECAClusters(0),fNTowerInGroup(0),fCalibData(0), + fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.), + fECAW0(0.),fRecPointsInRun(0),fTimeCut(0.),fMinECut(0.) { // default ctor (to be used mainly by Streamer) InitParameters() ; - fDefaultInit = kTRUE ; - for(Int_t is=0;is<20;is++){ - for(Int_t i=0;i<60;i++){ - for(Int_t j=0;j<60;j++){ - addOn[is][i][j]=0; - } - } - } -//PH cout<<"file to read 1"<GetTransformationForSM(); // Global <-> Local } //____________________________________________________________________________ AliEMCALClusterizerv1::AliEMCALClusterizerv1(const TString alirunFileName, const TString eventFolderName) -:AliEMCALClusterizer(alirunFileName, eventFolderName) + : AliEMCALClusterizer(alirunFileName, eventFolderName), + fHists(0),fPointE(0),fPointL1(0),fPointL2(0), + fPointDis(0),fPointMult(0),fDigitAmp(0),fMaxE(0), + fMaxL1(0),fMaxL2(0),fMaxDis(0),fGeom(0), + fDefaultInit(kFALSE), + fToUnfold(kFALSE), + fNumberOfECAClusters(0),fNTowerInGroup(0),fCalibData(0), + fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.), + fECAW0(0.),fRecPointsInRun(0),fTimeCut(0.),fMinECut(0.) { // ctor with the indication of the file where header Tree and digits Tree are stored InitParameters() ; Init() ; - fDefaultInit = kFALSE ; - for(Int_t is=0;is<20;is++){ - for(Int_t i=0;i<60;i++){ - for(Int_t j=0;j<60;j++){ - addOn[is][i][j]=0; - } - } - } -//PH cout<<"file to read 2"<GetCellIndex(AbsId, iSupMod, nModule, nIphi, nIeta) ; + if(!bCell) { + fGeom->PrintGeometry(); + Error("Calibrate()"," Wrong cell id number : %i", AbsId); + assert(0); + } + + fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta); + + fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi); + fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi); + + return -fADCpedestalECA + amp * fADCchannelECA ; + + } + else //Return energy with default parameters if calibration is not available + return -fADCpedestalECA + amp * fADCchannelECA ; + } //____________________________________________________________________________ void AliEMCALClusterizerv1::Exec(Option_t * option) { - // Steering method to perform clusterization for events - // in the range from fFirstEvent to fLastEvent. - // This range is optionally set by SetEventRange(). - // if fLastEvent=-1 (by default), then process events until the end. + // Steering method to perform clusterization for the current event + // in AliEMCALLoader if(strstr(option,"tim")) gBenchmark->Start("EMCALClusterizer"); if(strstr(option,"print")) Print("") ; + + AliRunLoader *rl = AliRunLoader::GetRunLoader(); + AliEMCALLoader *emcalLoader = dynamic_cast(rl->GetDetectorLoader("EMCAL")); - AliEMCALGetter * gime = AliEMCALGetter::Instance() ; + //Get calibration parameters from file or digitizer default values. + GetCalibrationParameters() ; - 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") ; - GetCalibrationParameters() ; + fNumberOfECAClusters = 0; - fNumberOfECAClusters = 0; - - MakeClusters() ; + if(strstr(option,"pseudo")) + MakeClusters("pseudo") ; //both types + else + MakeClusters("") ; //only the real clusters - if(fToUnfold) - MakeUnfolding() ; + if(fToUnfold) + MakeUnfolding() ; - WriteRecPoints() ; + WriteRecPoints() ; - if(strstr(option,"deb")) - PrintRecPoints(option) ; + if(strstr(option,"deb") || strstr(option,"all")) + PrintRecPoints(option) ; - //increment the total number of recpoints per run - fRecPointsInRun += gime->ECARecPoints()->GetEntriesFast() ; - } - - Unload(); + AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",emcalLoader->RecPoints()->GetEntriesFast())); + + //increment the total number of recpoints per run + fRecPointsInRun += emcalLoader->RecPoints()->GetEntriesFast() ; if(strstr(option,"tim")){ gBenchmark->Stop("EMCALClusterizer"); - printf("Exec took %f seconds for Clusterizing %f seconds per event", - gBenchmark->GetCpuTime("EMCALClusterizer"), gBenchmark->GetCpuTime("EMCALClusterizer")/nEvents ) ; - } - - SaveHists(); - + printf("Exec took %f seconds for Clusterizing", + gBenchmark->GetCpuTime("EMCALClusterizer")); + } } //____________________________________________________________________________ @@ -195,9 +258,9 @@ 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() ; - + 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) ; @@ -217,16 +280,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] ; @@ -281,16 +343,32 @@ 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(); -//PH cout<<"ChannelECA, peds "<GetStorage("local://CalibDB"); + + //Check if calibration is stored in data base + + AliEMCALLoader *emcalLoader = + dynamic_cast(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")); + + fCalibData =emcalLoader->CalibData(); + + 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(); + } } //____________________________________________________________________________ @@ -299,21 +377,19 @@ 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(); - if(!gime) - 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() ; -//PH cout<<"gime,geom "<GetTransformationForSM(); // Global <-> Local + AliInfo(Form("geom 0x%x",fGeom)); -//Sub fNTowers = geom->GetNZ() * geom->GetNPhi() ; - fNTowers =400; if(!gMinuit) gMinuit = new TMinuit(100) ; - //Sub if ( !gime->Clusterizer() ) - //Sub gime->PostClusterizer(this); - BookHists(); -//PH cout<<"hists booked "<EMCALGeometry() ; - - Int_t rv = 0 ; - - Int_t relid1[2] ; - //Sub geom->AbsToRelNumbering(d1->GetId(), relid1) ; - Int_t nSupMod=0; - Int_t nTower=0; - Int_t nIphi=0; - Int_t nIeta=0; - Int_t iphi=0; - Int_t ieta=0; - geom->GetCellIndex(d1->GetId(), nSupMod,nTower,nIphi,nIeta); - geom->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi,ieta); - relid1[0]=ieta; - relid1[1]=iphi; - - - Int_t nSupMod1=0; - Int_t nTower1=0; - Int_t nIphi1=0; - Int_t nIeta1=0; - Int_t iphi1=0; - Int_t ieta1=0; - geom->GetCellIndex(d2->GetId(), nSupMod1,nTower1,nIphi1,nIeta1); - geom->GetCellPhiEtaIndexInSModule(nSupMod1, nTower1,nIphi1,nIeta1, iphi1,ieta1); - Int_t relid2[2] ; - relid2[0]=ieta1; - relid2[1]=iphi1; - - //Sub geom->AbsToRelNumbering(d2->GetId(), relid2) ; - + 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 rowdiff = TMath::Abs( relid1[0] - relid2[0] ) ; - Int_t coldiff = TMath::Abs( relid1[1] - relid2[1] ) ; + fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1); + fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2); + if(nSupMod1 != nSupMod2) return 2; // different SM + + fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1); + fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2); + + rowdiff = TMath::Abs(iphi1 - iphi2); + coldiff = TMath::Abs(ieta1 - ieta2) ; - if (( coldiff <= 1 ) && ( rowdiff <= 1 )){ - rv = 1 ; - } - else { - if((relid2[0] > relid1[0]) && (relid2[1] > relid1[1]+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(coldiff + rowdiff <= 1) rv = 1; // neighbours with at least commom side; Nov 6,2006 - if (gDebug == 2 ) -if(rv==1)printf("AreNeighbours: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d \n", - 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() ; + // 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 + + //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 ; } //____________________________________________________________________________ @@ -410,163 +503,170 @@ 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(); + } + else if (treeR->GetEntries() > 0) { + Warning("WriteRecPoints","RecPoints already exist in output file. New Recpoitns will not be visible."); + } Int_t index ; //Evaluate position, dispersion and other RecPoint properties for EC section - for(index = 0; index < aECARecPoints->GetEntries(); index++) - (dynamic_cast(aECARecPoints->At(index)))->EvalAll(fECAW0,digits) ; + for(index = 0; index < aECARecPoints->GetEntries(); index++) { + if (dynamic_cast(aECARecPoints->At(index))->GetClusterType() != AliESDCaloCluster::kPseudoCluster) + dynamic_cast(aECARecPoints->At(index))->EvalAll(fECAW0,digits) ; + } aECARecPoints->Sort() ; - for(index = 0; index < aECARecPoints->GetEntries(); 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()); - - branchECA->Fill() ; + + TBranch * branchECA = 0; + if ((branchECA = treeR->GetBranch("EMCALECARP"))) + branchECA->SetAddress(&aECARecPoints); + else + treeR->Branch("EMCALECARP","TObjArray",&aECARecPoints,bufferSize,splitlevel); + treeR->Fill() ; - gime->WriteRecPoints("OVERWRITE"); - gime->WriteClusterizer("OVERWRITE"); + emcalLoader->WriteRecPoints("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 - - AliEMCALGetter * gime = AliEMCALGetter::Instance() ; - - AliEMCALGeometry * geom = gime->EMCALGeometry() ; + // Mar 03, 2007 by PAI - TObjArray * aECARecPoints = gime->ECARecPoints() ; + if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader"); - aECARecPoints->Delete() ; + AliEMCALLoader *emcalLoader = dynamic_cast(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")); + TObjArray * aECARecPoints = emcalLoader->RecPoints() ; + aECARecPoints->Clear(); - TClonesArray * digits = gime->Digits() ; - TClonesArray * digitsC = dynamic_cast(digits->Clone()) ; + TClonesArray *digits = emcalLoader->Digits(); - // Clusterization starts - TIter nextdigit(digitsC) ; - AliEMCALDigit * digit; + // Set up TObjArray with pointers to digits to work on + TObjArray *digitsC = new TObjArray(); + TIter nextdigit(digits); + AliEMCALDigit *digit; + while ( (digit = dynamic_cast(nextdigit())) ) { + digitsC->AddLast(digit); + } -//just for hist -/* - while ( (digit = dynamic_cast(nextdigit())) ) { // scan over the list of digitsC - digitAmp->Fill(digit->GetAmp()); + //Start with pseudoclusters, if option + if(strstr(option,"pseudo")) { + // + // New algorithm : will be created one pseudo cluster per module + // + AliDebug(1,Form("Pseudo clustering #digits : %i ",digits->GetEntries())); + + 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::kPseudoCluster); + } + recPoints[nSM]->AddDigit(*digit, Calibrate(digit->GetAmp(), digit->GetId())); + } + } + fNumberOfECAClusters = 0; + for(int i=0; iGetNumberOfSuperModules(); i++) { // put non empty rec.points to container + if(recPoints[i]) aECARecPoints->AddAt(recPoints[i], fNumberOfECAClusters++); + } + AliDebug(1,Form(" Number of PC %d ", fNumberOfECAClusters)); } - */ -/////////// - while ( (digit = dynamic_cast(nextdigit())) ) { // scan over the list of digitsC - AliEMCALRecPoint * clu = 0 ; - - TArrayI clusterECAdigitslist(50000); - -//Sub if (gDebug == 2) { - // printf("MakeClusters: id = %d, ene = %f , thre = %f \n", digit->GetId(),Calibrate(digit->GetAmp()), -// fECAClusteringThreshold) ; -// } -////////////////////////temp solution, embedding/////////////////// - int nSupMod=0, nTower=0, nIphi=0, nIeta=0; - int iphi=0, ieta=0; - geom->GetCellIndex(digit->GetId(), nSupMod,nTower,nIphi,nIeta); - geom->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi,ieta); - - /////////////////////////////////// - -//cout<IsInECA(digit->GetId()) && (Calibrate(digit->GetAmp()) > fECAClusteringThreshold ) ){ - if (geom->CheckAbsCellId(digit->GetId()) && (Calibrate(digit->GetAmp()+addOn[nSupMod-1][ieta-1][iphi-1]) > fECAClusteringThreshold ) ){ - //if(addOn[nSupMod-1][ieta-1][iphi-1]>0)cout<<"11 digit, add "<GetAmp()<(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", + digits->GetEntries(),fMinECut,ehs)); + + nextdigitC.Reset(); + + while ( (digit = dynamic_cast(nextdigitC())) ) { // scan over the list of digitsC + TArrayI clusterECAdigitslist(digits->GetEntries()); + + if(fGeom->CheckAbsCellId(digit->GetId()) && (Calibrate(digit->GetAmp(), digit->GetId()) > fECAClusteringThreshold ) ){ // 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)) ; + if(fNumberOfECAClusters >= aECARecPoints->GetSize()) aECARecPoints->Expand(2*fNumberOfECAClusters+1) ; + AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("") ; + aECARecPoints->AddAt(recPoint, fNumberOfECAClusters) ; + recPoint = dynamic_cast(aECARecPoints->At(fNumberOfECAClusters)) ; fNumberOfECAClusters++ ; - clu->AddDigit(*digit, Calibrate(digit->GetAmp()+addOn[nSupMod-1][ieta-1][iphi-1])) ; - clusterECAdigitslist[iDigitInECACluster] = digit->GetIndexInList() ; - iDigitInECACluster++ ; -// cout<<" 1st setp:cluno, digNo "<SetClusterType(AliESDCaloCluster::kClusterv1); + + recPoint->AddDigit(*digit, Calibrate(digit->GetAmp(), digit->GetId())) ; + TObjArray clusterDigits; + clusterDigits.AddLast(digit); digitsC->Remove(digit) ; - if (gDebug == 2 ) - printf("MakeClusters: OK id = %d, ene = %f , thre = %f \n", digit->GetId(),Calibrate(digit->GetAmp()), fECAClusteringThreshold) ; - nextdigit.Reset() ; + + 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 ; - - // 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 -// cout<<"we have new digit "<GetCellIndex(digitN->GetId(), nSupMod,nTower,nIphi,nIeta); - geom->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, iphi,ieta); - //////////////// - if(Calibrate(digitN->GetAmp()+addOn[nSupMod-1][ieta-1][iphi-1]) < fMinECut ) digitsC->Remove(digitN); -// cout<<" new digit above ECut "<0)cout<<"22 digit, add "<GetAmp()<AddDigit(*digitN, Calibrate( digitN->GetAmp()+addOn[nSupMod-1][ieta-1][iphi-1]) ) ; - clusterECAdigitslist[iDigitInECACluster] = digitN->GetIndexInList() ; - iDigitInECACluster++ ; -// cout<<"when neighbour: cluno, digNo "<GetId()<<" "<Remove(digitN) ; -// break ; -// case 2 : // too far from each other -//Subh goto endofloop1; -// cout<<"earlier go to end of loop"<GetAmp()+addOn[nSupMod-1][ieta-1][iphi-1]) < fMinECut ){ - //if(addOn[nSupMod-1][ieta-1][iphi-1]>0)cout<<"33 digit, add "<GetAmp()<Remove(digit); - } - //cout<<"after endofloop: cluno, digNo "<(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 ; -cout<<"total no of clusters "<GetEntriesFast()<<" digits"<GetEntriesFast())); } -//____________________________________________________________________________ void AliEMCALClusterizerv1::MakeUnfolding() const { Fatal("AliEMCALClusterizerv1::MakeUnfolding", "--> Unfolding not implemented") ; - } //____________________________________________________________________________ @@ -643,26 +743,33 @@ 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:") ; + AliEMCALLoader *emcalLoader = dynamic_cast(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")); + TObjArray * aECARecPoints = emcalLoader->RecPoints() ; + + if(strstr(option,"deb")) { + printf("PrintRecPoints: Clusterization result:") ; - printf("event # %d\n", gAlice->GetEvNumber() ) ; - printf(" Found %d ECA Rec Points\n ", + printf("event # %d\n", emcalLoader->GetRunLoader()->GetEventNumber() ) ; + printf(" Found %d ECA Rec Points\n ", aECARecPoints->GetEntriesFast()) ; + } fRecPointsInRun += aECARecPoints->GetEntriesFast() ; if(strstr(option,"all")) { - Int_t index =0; - printf("\n-----------------------------------------------------------------------\n") ; - printf("Clusters in ECAL section\n") ; - printf("Index Ene(GeV) Multi Module phi r theta X Y Z Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ; + if(strstr(option,"deb")) { + printf("\n-----------------------------------------------------------------------\n") ; + printf("Clusters in ECAL section\n") ; + printf("Index Ene(GeV) Multi Module GX GY GZ lX lY lZ Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ; + } + Int_t index =0; Float_t maxE=0; Float_t maxL1=0; Float_t maxL2=0; Float_t maxDis=0; + AliEMCALHistoUtilities::FillH1(fHists, 12, double(aECARecPoints->GetEntries())); + for (index = 0 ; index < aECARecPoints->GetEntries() ; index++) { AliEMCALRecPoint * rp = dynamic_cast(aECARecPoints->At(index)) ; TVector3 globalpos; @@ -674,6 +781,7 @@ void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option) Int_t * primaries; Int_t nprimaries; primaries = rp->GetPrimaries(nprimaries); + if(strstr(option,"deb")) printf("\n%6d %8.4f %3d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ", rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(), globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(), @@ -685,66 +793,85 @@ void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option) maxL2=lambda[1]; maxDis=rp->GetDispersion(); } - pointE->Fill(rp->GetEnergy()); - pointL1->Fill(lambda[0]); - pointL2->Fill(lambda[1]); - pointDis->Fill(rp->GetDispersion()); - pointMult->Fill(rp->GetMultiplicity()); + 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); - MaxL1->Fill(maxL1); - MaxL2->Fill(maxL2); - MaxDis->Fill(maxDis); - + fMaxE->Fill(maxE); + fMaxL1->Fill(maxL1); + fMaxL2->Fill(maxL2); + fMaxDis->Fill(maxDis); + if(strstr(option,"deb")) printf("\n-----------------------------------------------------------------------\n"); } } - void AliEMCALClusterizerv1::BookHists() +TList* AliEMCALClusterizerv1::BookHists() { - pointE = new TH1F("pointE","point energy", 2000, 0.0, 150.); - pointL1 = new TH1F("pointL1","point L1", 1000, 0.0, 3.); - pointL2 = new TH1F("pointL2","point L2", 1000, 0.0, 3.); - pointDis = new TH1F("pointDis","point Dis", 1000, 0.0, 3.); - pointMult = new TH1F("pointMult","point Mult", 100, 0.0, 100.); - digitAmp = new TH1F("digitAmp","Digit Amplitude", 2000, 0.0, 5000.); - MaxE = new TH1F("maxE","Max point energy", 2000, 0.0, 150.); - MaxL1 = new TH1F("maxL1","Max point L1", 1000, 0.0, 3.); - MaxL2 = new TH1F("maxL2","Max point L2", 1000, 0.0, 3.); - MaxDis = new TH1F("maxDis","Max point Dis", 1000, 0.0, 3.); + //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() + +void AliEMCALClusterizerv1::SaveHists(const char *fn) { - recofile=new TFile("reco.root","RECREATE"); - pointE->Write(); - pointL1->Write(); - pointL2->Write(); - pointDis->Write(); - pointMult->Write(); - digitAmp->Write(); - MaxE->Write(); - MaxL1->Write(); - MaxL2->Write(); - MaxDis->Write(); -recofile->Close(); + AliEMCALHistoUtilities::SaveListOfHists(fHists, fn, kTRUE); } -void AliEMCALClusterizerv1::ReadFile() +void AliEMCALClusterizerv1::PrintRecoInfo() { - return; // 3-jan-05 - FILE *fp = fopen("hijing1.dat","r"); - for(Int_t line=0;line<9113;line++){ - Int_t eg,l1,l2,sm; - Int_t ncols0; - ncols0 = fscanf(fp,"%d %d %d %d",&sm,&l1,&l2,&eg); - // cout<At(12); + if(h) { + printf(" ## Multiplicity of RecPoints ## \n"); + for(int i=1; i<=h->GetNbinsX(); i++) { + int nbin = int((*h)[i]); + int mult = int(h->GetBinCenter(i)); + if(nbin > 0) printf(" %i : %5.5i %6.3f %% \n", mult, nbin, 100.*nbin/h->GetEntries()); + } + } } +void AliEMCALClusterizerv1::DrawLambdasHists() +{ + if(fMaxL1) { + fMaxL1->Draw(); + if(fMaxL2) fMaxL2->Draw("same"); + if(fMaxDis) { + fMaxDis->Draw("same"); + } + } +} + +void AliEMCALClusterizerv1::Browse(TBrowser* b) +{ + if(fHists) b->Add(fHists); + if(fGeom) b->Add(fGeom); + TTask::Browse(b); +}