X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliEMCALClusterizerv1.cxx;h=a654112c454dd6dad4abc5e66d4e39f2a470024f;hb=b85ea106fc895e812397f1f3bab3ce9d20798225;hp=a298ffefbcf4637037d7f4722407fb025323c8cb;hpb=2bb3725c6f8db942aa39400603f72c38f6d97be0;p=u%2Fmrichter%2FAliRoot.git diff --git a/EMCAL/AliEMCALClusterizerv1.cxx b/EMCAL/AliEMCALClusterizerv1.cxx index a298ffefbcf..a654112c454 100644 --- a/EMCAL/AliEMCALClusterizerv1.cxx +++ b/EMCAL/AliEMCALClusterizerv1.cxx @@ -15,9 +15,10 @@ /* $Id$ */ -//*-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute) +//-- Author: Yves Schutz (SUBATECH) & Dmitri Peressounko (SUBATECH & Kurchatov Institute) // August 2002 Yves Schutz: clone PHOS as closely as possible and intoduction // of new IO (à la PHOS) +// Mar 2007, Aleksei Pavlinov - new algoritmh of pseudo clusters ////////////////////////////////////////////////////////////////////////////// // Clusterization class. Performs clusterization (collects neighbouring active cells) and // unfolds the clusters having several local maxima. @@ -45,6 +46,7 @@ // // time - print benchmarking results // --- ROOT system --- +#include class TROOT; #include @@ -56,6 +58,7 @@ class TFolder; class TSystem; #include #include +#include // --- Standard library --- @@ -64,14 +67,14 @@ class TSystem; #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 "AliEMCALHistoUtilities.h" +#include "AliEMCALRecParam.h" +#include "AliEMCALReconstructor.h" #include "AliCDBManager.h" class AliCDBStorage; @@ -80,86 +83,76 @@ class AliCDBStorage; ClassImp(AliEMCALClusterizerv1) //____________________________________________________________________________ -AliEMCALClusterizerv1::AliEMCALClusterizerv1() +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), + fGeom(0), + fDefaultInit(kFALSE), fToUnfold(kFALSE), - fNumberOfECAClusters(0),fNTowerInGroup(0),fCalibData(0), + fNumberOfECAClusters(0),fCalibData(0), fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.), - fECAW0(0.),fRecPointsInRun(0),fTimeGate(0.),fMinECut(0.) + fECAW0(0.),fTimeCut(0.),fMinECut(0.) { - // default ctor (to be used mainly by Streamer) + // ctor with the indication of the file where header Tree and digits Tree are stored - InitParameters() ; - fGeom = AliEMCALGeometry::GetInstance(); - fGeom->GetTransformationForSM(); // Global <-> Local + Init() ; } //____________________________________________________________________________ -AliEMCALClusterizerv1::AliEMCALClusterizerv1(const TString alirunFileName, const TString 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), +AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry) + : AliEMCALClusterizer(), + fGeom(geometry), fDefaultInit(kFALSE), fToUnfold(kFALSE), - fNumberOfECAClusters(0),fNTowerInGroup(0),fCalibData(0), + fNumberOfECAClusters(0),fCalibData(0), fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.), - fECAW0(0.),fRecPointsInRun(0),fTimeGate(0.),fMinECut(0.) + fECAW0(0.),fTimeCut(0.),fMinECut(0.) { // ctor with the indication of the file where header Tree and digits Tree are stored - - InitParameters() ; - Init() ; + // use this contructor to avoid usage of Init() which uses runloader + // change needed by HLT - MP + + // Note for the future: the use on runloader should be avoided or optional at least + // another way is to make Init virtual and protected at least such that the deriving classes can overload + // Init() ; + // + + if (!fGeom) + { + AliFatal("Geometry not initialized."); + } + + if(!gMinuit) + gMinuit = new TMinuit(100) ; + } //____________________________________________________________________________ -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) +AliEMCALClusterizerv1::AliEMCALClusterizerv1(AliEMCALGeometry* geometry, AliEMCALCalibData * calib) +: AliEMCALClusterizer(), +fGeom(geometry), +fDefaultInit(kFALSE), +fToUnfold(kFALSE), +fNumberOfECAClusters(0),fCalibData(calib), +fADCchannelECA(0.),fADCpedestalECA(0.),fECAClusteringThreshold(0.),fECALocMaxCut(0.), +fECAW0(0.),fTimeCut(0.),fMinECut(0.) { - //copy ctor + // ctor, geometry and calibration are initialized elsewhere. + + if (!fGeom) + AliFatal("Geometry not initialized."); + + if(!gMinuit) + gMinuit = new TMinuit(100) ; + } + //____________________________________________________________________________ AliEMCALClusterizerv1::~AliEMCALClusterizerv1() { // dtor } -//____________________________________________________________________________ -const TString AliEMCALClusterizerv1::BranchName() const -{ - return GetName(); -} - //____________________________________________________________________________ Float_t AliEMCALClusterizerv1::Calibrate(Int_t amp, Int_t AbsId) { @@ -169,18 +162,6 @@ Float_t AliEMCALClusterizerv1::Calibrate(Int_t amp, Int_t AbsId) // or from digitizer parameters for simulated data. if(fCalibData){ - - //JLK 13-Mar-2006 - //We now get geometry at a higher level - // - // Loader - // AliRunLoader *rl = AliRunLoader::GetRunLoader(); - - // Load EMCAL Geomtry - // rl->LoadgAlice(); - //AliRun * gAlice = rl->GetAliRun(); - //AliEMCAL * emcal = (AliEMCAL*)gAlice->GetDetector("EMCAL"); - //AliEMCALGeometry * geom = emcal->GetGeometry(); if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader") ; @@ -198,11 +179,13 @@ Float_t AliEMCALClusterizerv1::Calibrate(Int_t amp, Int_t AbsId) Error("Calibrate()"," Wrong cell id number : %i", AbsId); assert(0); } + fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta); fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi); fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi); - return -fADCpedestalECA + amp * fADCchannelECA ; + + return -fADCpedestalECA + amp * fADCchannelECA ; } else //Return energy with default parameters if calibration is not available @@ -211,147 +194,140 @@ Float_t AliEMCALClusterizerv1::Calibrate(Int_t amp, Int_t AbsId) } //____________________________________________________________________________ -void AliEMCALClusterizerv1::Exec(Option_t * option) +void AliEMCALClusterizerv1::Digits2Clusters(Option_t * option) { - // Steering method to perform clusterization for events - // in the range from fFirstEvent to fLastEvent. - // This range is optionally set by SetEventRange(). - // if fLastEvent=-1 (by default), then process events until the end. + // Steering method to perform clusterization for the current event + // in AliEMCALLoader if(strstr(option,"tim")) gBenchmark->Start("EMCALClusterizer"); if(strstr(option,"print")) Print("") ; - - AliRunLoader *rl = AliRunLoader::GetRunLoader(); - AliEMCALLoader *emcalLoader = dynamic_cast(rl->GetDetectorLoader("EMCAL")); - + //Get calibration parameters from file or digitizer default values. GetCalibrationParameters() ; - if (fLastEvent == -1) - fLastEvent = rl->GetNumberOfEvents() - 1; - Int_t nEvents = fLastEvent - fFirstEvent + 1; - Int_t ievent ; - rl->LoadDigits("EMCAL"); - for (ievent = fFirstEvent; ievent <= fLastEvent; ievent++) { - rl->GetEvent(ievent); + fNumberOfECAClusters = 0; - fNumberOfECAClusters = 0; + MakeClusters() ; //only the real clusters - if(strstr(option,"pseudo")) - MakeClusters("pseudo") ; //both types - else - MakeClusters("") ; //only the real clusters + if(fToUnfold) + MakeUnfolding() ; - if(fToUnfold) - MakeUnfolding() ; + Int_t index ; - WriteRecPoints() ; + //Evaluate position, dispersion and other RecPoint properties for EC section + for(index = 0; index < fRecPoints->GetEntries(); index++) { + dynamic_cast(fRecPoints->At(index))->EvalAll(fECAW0,fDigitsArr) ; + } - if(strstr(option,"deb") || strstr(option,"all")) - PrintRecPoints(option) ; + fRecPoints->Sort() ; - //increment the total number of recpoints per run - fRecPointsInRun += emcalLoader->RecPoints()->GetEntriesFast() ; + for(index = 0; index < fRecPoints->GetEntries(); index++) { + (dynamic_cast(fRecPoints->At(index)))->SetIndexInList(index) ; + (dynamic_cast(fRecPoints->At(index)))->Print(); } + + fTreeR->Fill(); - Unload(); + if(strstr(option,"deb") || strstr(option,"all")) + PrintRecPoints(option) ; + + AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast())); if(strstr(option,"tim")){ gBenchmark->Stop("EMCALClusterizer"); - printf("Exec took %f seconds for Clusterizing %f seconds per event", - gBenchmark->GetCpuTime("EMCALClusterizer"), gBenchmark->GetCpuTime("EMCALClusterizer")/nEvents ); - } - + printf("Exec took %f seconds for Clusterizing", + gBenchmark->GetCpuTime("EMCALClusterizer")); + } } //____________________________________________________________________________ -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 +Bool_t AliEMCALClusterizerv1::FindFit(AliEMCALRecPoint * RecPoint, 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 - AliEMCALLoader *emcalLoader = dynamic_cast(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")); - TClonesArray *digits = emcalLoader->Digits(); + if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader"); gMinuit->mncler(); // Reset Minuit's list of paramters gMinuit->SetPrintLevel(-1) ; // No Printout - gMinuit->SetFCN(AliEMCALClusterizerv1::UnfoldingChiSquare) ; - // To set the address of the minimization function + gMinuit->SetFCN(AliEMCALClusterizerv1::UnfoldingChiSquare) ; + // To set the address of the minimization function TList * toMinuit = new TList(); - toMinuit->AddAt(emcRP,0) ; - toMinuit->AddAt(digits,1) ; - + toMinuit->AddAt(RecPoint,0) ; + toMinuit->AddAt(fDigitsArr,1) ; + toMinuit->AddAt(fGeom,2) ; + gMinuit->SetObjectFit(toMinuit) ; // To tranfer pointer to UnfoldingChiSquare // filling initial values for fit parameters AliEMCALDigit * digit ; - Int_t ierflg = 0; + Int_t ierflg = 0; Int_t index = 0 ; Int_t nDigits = (Int_t) nPar / 3 ; Int_t iDigit ; for(iDigit = 0; iDigit < nDigits; iDigit++){ - digit = maxAt[iDigit]; + digit = maxAt[iDigit]; + Double_t x = 0.; + Double_t y = 0.; + Double_t z = 0.; - Float_t x = 0.; - Float_t z = 0.; - // 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 + fGeom->RelPosCellInSModule(digit->GetId(), y, x, z); Float_t energy = maxAtEnergy[iDigit] ; gMinuit->mnparm(index, "x", x, 0.1, 0, 0, ierflg) ; - index++ ; - if(ierflg != 0){ - Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : x = %f", x ) ; + index++ ; + if(ierflg != 0){ + Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : x = %f", x ) ; return kFALSE; } gMinuit->mnparm(index, "z", z, 0.1, 0, 0, ierflg) ; - index++ ; + index++ ; if(ierflg != 0){ - Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : z = %f", z) ; + Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : z = %f", z) ; return kFALSE; } gMinuit->mnparm(index, "Energy", energy , 0.05*energy, 0., 4.*energy, ierflg) ; - index++ ; + index++ ; if(ierflg != 0){ - Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : energy = %f", energy) ; + Error("FindFit", "EMCAL Unfolding Unable to set initial value for fit procedure : energy = %f", energy) ; return kFALSE; } } - Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; The number of function call slightly - // depends on it. - Double_t p1 = 1.0 ; + Double_t p0 = 0.1 ; // "Tolerance" Evaluation stops when EDM = 0.0001*p0 ; + // The number of function call slightly depends on it. + //Double_t p1 = 1.0 ; Double_t p2 = 0.0 ; - gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls - gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient + gMinuit->mnexcm("SET STR", &p2, 0, ierflg) ; // force TMinuit to reduce function calls + // gMinuit->mnexcm("SET GRA", &p1, 1, ierflg) ; // force TMinuit to use my gradient gMinuit->SetMaxIterations(5); gMinuit->mnexcm("SET NOW", &p2 , 0, ierflg) ; // No Warnings + gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize - gMinuit->mnexcm("MIGRAD", &p0, 0, ierflg) ; // minimize - - if(ierflg == 4){ // Minimum not found - Error("FindFit", "EMCAL Unfolding Fit not converged, cluster abandoned " ) ; + if(ierflg == 4){ // Minimum not found + Error("FindFit", "EMCAL Unfolding Fit not converged, cluster abandoned " ) ; return kFALSE ; - } + } for(index = 0; index < nPar; index++){ Double_t err ; Double_t val ; gMinuit->GetParameter(index, val, err) ; // Returns value and error of parameter index fitparameters[index] = val ; - } + } delete toMinuit ; return kTRUE; @@ -370,24 +346,17 @@ void AliEMCALClusterizerv1::GetCalibrationParameters() // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB"); //Check if calibration is stored in data base - if(AliCDBManager::Instance()->IsDefaultStorageSet()){ - AliCDBEntry *entry = (AliCDBEntry*) - AliCDBManager::Instance()->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(); - } + + if(!fCalibData) + { + AliCDBEntry *entry = (AliCDBEntry*) + AliCDBManager::Instance()->Get("EMCAL/Calib/Data"); + if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject(); + } + + if(!fCalibData) + AliFatal("Calibration parameters not found in CDB!"); + } //____________________________________________________________________________ @@ -396,15 +365,17 @@ 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 - AliRunLoader *rl = AliRunLoader::GetRunLoader(); - fGeom = dynamic_cast(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry(); - fGeom->GetTransformationForSM(); // Global <-> Local - AliInfo(Form("geom 0x%x",fGeom)); + AliRunLoader *rl = AliRunLoader::Instance(); + if (rl->GetAliRun() && rl->GetAliRun()->GetDetector("EMCAL")) + fGeom = dynamic_cast(rl->GetAliRun()->GetDetector("EMCAL"))->GetGeometry(); + else + fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName()); + + AliDebug(1,Form("geom 0x%x",fGeom)); if(!gMinuit) gMinuit = new TMinuit(100) ; - fHists = BookHists(); } //____________________________________________________________________________ @@ -412,19 +383,25 @@ void AliEMCALClusterizerv1::InitParameters() { // Initializes the parameters for the Clusterizer fNumberOfECAClusters = 0; + fTimeCut = 300e-9 ; // 300 ns time cut (to be tuned) - fNTowerInGroup = 36; //Produces maximum of 80 pseudoclusters per event - - fECAClusteringThreshold = 0.1; // value obtained from Aleksei - fECALocMaxCut = 0.03; // ?? + fCalibData = 0 ; - fECAW0 = 4.5; - fTimeGate = 1.e-8 ; - fToUnfold = kFALSE ; - fRecPointsInRun = 0 ; - fMinECut = 0.01; // have to be tune + const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam(); + if(!recParam) { + AliFatal("Reconstruction parameters for EMCAL not set!"); + } else { + fECAClusteringThreshold = recParam->GetClusteringThreshold(); + fECAW0 = recParam->GetW0(); + fMinECut = recParam->GetMinECut(); + fToUnfold = recParam->GetUnfold(); + if(fToUnfold) AliWarning("Cluster Unfolding ON. Implementing only for eta=0 case!!!"); + fECALocMaxCut = recParam->GetLocMaxCut(); + + AliDebug(1,Form("Reconstruction parameters: fECAClusteringThreshold=%.3f, fECAW=%.3f, fMinECut=%.3f, fToUnfold=%d, fECALocMaxCut=%.3f", + fECAClusteringThreshold,fECAW0,fMinECut,fToUnfold,fECALocMaxCut)); + } - fCalibData = 0 ; } //____________________________________________________________________________ @@ -453,8 +430,8 @@ Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d rowdiff = TMath::Abs(iphi1 - iphi2); coldiff = TMath::Abs(ieta1 - ieta2) ; - // if (( coldiff <= 1 ) && ( rowdiff <= 1 )) rv = 1; // neighbours with at least commom vertex - if(coldiff + rowdiff <= 1) rv = 1; // neighbours with at least commom side; Nov 6,2006 + // neighbours with at least commom side; May 11, 2007 + if ((coldiff==0 && abs(rowdiff)==1) || (rowdiff==0 && abs(coldiff)==1)) rv = 1; if (gDebug == 2 && rv==1) printf("AreNeighbours: neighbours=%d, id1=%d, relid1=%d,%d \n id2=%d, relid2=%d,%d \n", @@ -464,318 +441,329 @@ Int_t AliEMCALClusterizerv1::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d } //____________________________________________________________________________ -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, 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 ; -} - -//____________________________________________________________________________ -void AliEMCALClusterizerv1::Unload() -{ - // Unloads the Digits and RecPoints - AliEMCALLoader *emcalLoader = dynamic_cast(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")); - - emcalLoader->UnloadDigits() ; - emcalLoader->UnloadRecPoints() ; -} - -//____________________________________________________________________________ -void AliEMCALClusterizerv1::WriteRecPoints() -{ - - // Creates new branches with given title - // fills and writes into TreeR. - - AliEMCALLoader *emcalLoader = dynamic_cast(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")); - - TObjArray * aECARecPoints = emcalLoader->RecPoints() ; - - 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) ; - - aECARecPoints->Sort() ; - - for(index = 0; index < aECARecPoints->GetEntries(); index++) { - (dynamic_cast(aECARecPoints->At(index)))->SetIndexInList(index) ; - (dynamic_cast(aECARecPoints->At(index)))->Print(); - } - - Int_t bufferSize = 32000 ; - Int_t splitlevel = 0 ; - - //EC section branch - - TBranch * branchECA = 0; - if ((branchECA = treeR->GetBranch("EMCALECARP"))) - branchECA->SetAddress(&aECARecPoints); - else - treeR->Branch("EMCALECARP","TObjArray",&aECARecPoints,bufferSize,splitlevel); - treeR->Fill() ; - - emcalLoader->WriteRecPoints("OVERWRITE"); - -} - -//____________________________________________________________________________ -void AliEMCALClusterizerv1::MakeClusters(char* option) +void AliEMCALClusterizerv1::MakeClusters() { // Steering method to construct the clusters stored in a list of Reconstructed Points // A cluster is defined as a list of neighbour digits + // Mar 03, 2007 by PAI - AliEMCALLoader *emcalLoader = dynamic_cast(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")); - - TObjArray * aECARecPoints = emcalLoader->RecPoints() ; - - if (fGeom==0) - AliFatal("Did not get geometry from EMCALLoader"); - - aECARecPoints->Clear(); - - //Start with pseudoclusters, if option - if(strstr(option,"pseudo")) { // ?? Nov 3, 2006 - TClonesArray * digits = emcalLoader->Digits() ; - TClonesArray * digitsC = dynamic_cast(digits->Clone()); - - TIter nextdigit(digitsC) ; - AliEMCALDigit * digit; - - AliEMCALRecPoint * recPoint = 0 ; - int ineb=0; - nextdigit.Reset(); - - // PseudoClusterization starts - while ( (digit = dynamic_cast(nextdigit())) ) { // scan over the list of digitsC - recPoint = 0 ; - TArrayI clusterECAdigitslist(1000); // what is this - - if(fGeom->CheckAbsCellId(digit->GetId()) ) { //Is this an EMCAL digit? Just maing sure... - Int_t iDigitInECACluster = 0; // start a new Tower RecPoint - if(fNumberOfECAClusters >= aECARecPoints->GetSize()) aECARecPoints->Expand(2*fNumberOfECAClusters+1) ; - AliEMCALRecPoint * rp = new AliEMCALRecPoint("") ; - aECARecPoints->AddAt(rp, fNumberOfECAClusters) ; - recPoint = dynamic_cast(aECARecPoints->At(fNumberOfECAClusters)) ; - fNumberOfECAClusters++ ; - - recPoint->SetClusterType(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(recPoint) - AliDebug(1,Form("MakeClusters: cl.e %f \n", recPoint->GetEnergy())); - //if(recPoint) cout << "cl.e " << recPoint->GetEnergy() << endl; - delete digitsC ; - } - - //Now do real clusters - TClonesArray * digits = emcalLoader->Digits() ; - TClonesArray * digitsC = dynamic_cast(digits->Clone()); // will work with thic copy + if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader"); - TIter nextdigitC(digitsC) ; - AliEMCALDigit * digit; + fRecPoints->Clear(); - AliEMCALRecPoint * recPoint = 0 ; - int ineb=0; - nextdigitC.Reset(); + // Set up TObjArray with pointers to digits to work on + TObjArray *digitsC = new TObjArray(); + TIter nextdigit(fDigitsArr); + AliEMCALDigit *digit; + while ( (digit = dynamic_cast(nextdigit())) ) { + digitsC->AddLast(digit); + } double e = 0.0, ehs = 0.0; + TIter nextdigitC(digitsC); + while ( (digit = dynamic_cast(nextdigitC())) ) { // clean up digits e = Calibrate(digit->GetAmp(), digit->GetId()); - AliEMCALHistoUtilities::FillH1(fHists, 10, digit->GetAmp()); - AliEMCALHistoUtilities::FillH1(fHists, 11, e); - if(e < fMinECut ) digitsC->Remove(digit); - else ehs += 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)); - //cout << " Number of digits " << digits->GetEntries() << " -> (e>" <GetEntries()<< " ehs "<GetEntries(),fMinECut,ehs)); - // Clusterization starts - // cout << "Outer Loop" << endl; - ineb=0; nextdigitC.Reset(); - recPoint = 0 ; + while ( (digit = dynamic_cast(nextdigitC())) ) { // scan over the list of digitsC - recPoint = 0 ; - TArrayI clusterECAdigitslist(1000); // what is this + TArrayI clusterECAdigitslist(fDigitsArr->GetEntries()); 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)) ; + // start a new Tower RecPoint + if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(2*fNumberOfECAClusters+1) ; + + AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("") ; + fRecPoints->AddAt(recPoint, fNumberOfECAClusters) ; + recPoint = dynamic_cast(fRecPoints->At(fNumberOfECAClusters)) ; fNumberOfECAClusters++ ; - recPoint->SetClusterType(AliESDCaloCluster::kClusterv1); + recPoint->SetClusterType(AliESDCaloCluster::kEMCALClusterv1); recPoint->AddDigit(*digit, Calibrate(digit->GetAmp(), digit->GetId())) ; - clusterECAdigitslist[iDigitInECACluster] = digit->GetIndexInList() ; - iDigitInECACluster++ ; + TObjArray clusterDigits; + clusterDigits.AddLast(digit); digitsC->Remove(digit) ; - AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , thre = %f \n", digit->GetId(), + + AliDebug(1,Form("MakeClusters: OK id = %d, ene = %f , cell.th. = %f \n", digit->GetId(), Calibrate(digit->GetAmp(),digit->GetId()), fECAClusteringThreshold)); - nextdigitC.Reset(); // will start from beggining - AliEMCALDigit * digitN = 0; // digi neighbor - Int_t index = 0 ; - - // Find the neighbours - NEIGHBOURS: - while (index < iDigitInECACluster){ // scan over digits already in cluster - digit = (AliEMCALDigit*)digits->At(clusterECAdigitslist[index]); - index++ ; - while ( (digitN = (AliEMCALDigit *)nextdigitC())) { // scan over the reduced list of digits - - ineb = AreNeighbours(digit, digitN); // call (digit,digitN) in THAT oder !!!!! - if(ineb==1) { + // Grow cluster by finding neighbours + TIter nextClusterDigit(&clusterDigits); + while ( (digit = dynamic_cast(nextClusterDigit())) ) { // scan over digits in cluster + TIter nextdigitN(digitsC); + AliEMCALDigit *digitN = 0; // digi neighbor + while ( (digitN = (AliEMCALDigit *)nextdigitN()) ) { // scan over all digits to look for neighbours + if (AreNeighbours(digit, digitN)==1) { // call (digit,digitN) in THAT oder !!!!! recPoint->AddDigit(*digitN, Calibrate(digitN->GetAmp(),digitN->GetId()) ) ; - clusterECAdigitslist[iDigitInECACluster] = digitN->GetIndexInList() ; - iDigitInECACluster++ ; + clusterDigits.AddLast(digitN) ; digitsC->Remove(digitN) ; - // Have to start from begining for clusters digits too ; Nov 3,2006 - index = 0; - nextdigitC.Reset() ; - goto NEIGHBOURS; } // if(ineb==1) - - } // scan over the reduced list of digits - nextdigitC.Reset() ; // will start from beginning + } // 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 - 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())); + delete digitsC ; + + AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast())); } //____________________________________________________________________________ -void AliEMCALClusterizerv1::MakeUnfolding() const +void AliEMCALClusterizerv1::MakeUnfolding() { - Fatal("AliEMCALClusterizerv1::MakeUnfolding", "--> Unfolding not implemented") ; - + // Unfolds clusters using the shape of an ElectroMagnetic shower + // Performs unfolding of all clusters + + if(fNumberOfECAClusters > 0){ + if (fGeom==0) + AliFatal("Did not get geometry from EMCALLoader") ; + Int_t nModulesToUnfold = fGeom->GetNCells(); + + Int_t numberofNotUnfolded = fNumberOfECAClusters ; + Int_t index ; + for(index = 0 ; index < numberofNotUnfolded ; index++){ + + AliEMCALRecPoint * RecPoint = dynamic_cast( fRecPoints->At(index) ) ; + + TVector3 gpos; + Int_t absId; + RecPoint->GetGlobalPosition(gpos); + fGeom->GetAbsCellIdFromEtaPhi(gpos.Eta(),gpos.Phi(),absId); + if(absId > nModulesToUnfold) + break ; + + Int_t nMultipl = RecPoint->GetMultiplicity() ; + AliEMCALDigit ** maxAt = new AliEMCALDigit*[nMultipl] ; + Float_t * maxAtEnergy = new Float_t[nMultipl] ; + Int_t nMax = RecPoint->GetNumberOfLocalMax(maxAt, maxAtEnergy,fECALocMaxCut,fDigitsArr) ; + + if( nMax > 1 ) { // if cluster is very flat (no pronounced maximum) then nMax = 0 + UnfoldCluster(RecPoint, nMax, maxAt, maxAtEnergy) ; + fRecPoints->Remove(RecPoint); + fRecPoints->Compress() ; + index-- ; + fNumberOfECAClusters-- ; + numberofNotUnfolded-- ; + } + else{ + RecPoint->SetNExMax(1) ; //Only one local maximum + } + + delete[] maxAt ; + delete[] maxAtEnergy ; + } + } + // End of Unfolding of clusters } //____________________________________________________________________________ -Double_t AliEMCALClusterizerv1::ShowerShape(Double_t r) +Double_t AliEMCALClusterizerv1::ShowerShape(Double_t x, Double_t y) { - // Shape of the shower (see EMCAL TDR) + // Shape of the shower // If you change this function, change also the gradient evaluation in ChiSquare() - Double_t r4 = r*r*r*r ; - Double_t r295 = TMath::Power(r, 2.95) ; - Double_t shape = TMath::Exp( -r4 * (1. / (2.32 + 0.26 * r4) + 0.0316 / (1 + 0.0652 * r295) ) ) ; + Double_t r = sqrt(x*x+y*y); + Double_t r133 = TMath::Power(r, 1.33) ; + Double_t r669 = TMath::Power(r, 6.69) ; + Double_t shape = TMath::Exp( -r133 * (1. / (1.57 + 0.0860 * r133) - 0.55 / (1 + 0.000563 * r669) ) ) ; return shape ; } //____________________________________________________________________________ -void AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * /*iniTower*/, - Int_t /*nMax*/, - AliEMCALDigit ** /*maxAt*/, - Float_t * /*maxAtEnergy*/) const +void AliEMCALClusterizerv1::UnfoldCluster(AliEMCALRecPoint * iniTower, + Int_t nMax, + AliEMCALDigit ** maxAt, + Float_t * maxAtEnergy) { // Performs the unfolding of a cluster with nMax overlapping showers - - Fatal("UnfoldCluster", "--> Unfolding not implemented") ; + Int_t nPar = 3 * nMax ; + Float_t * fitparameters = new Float_t[nPar] ; + + if (fGeom==0) + AliFatal("Did not get geometry from EMCALLoader") ; + + Bool_t rv = FindFit(iniTower, maxAt, maxAtEnergy, nPar, fitparameters) ; + if( !rv ) { + // Fit failed, return and remove cluster + iniTower->SetNExMax(-1) ; + delete[] fitparameters ; + return ; + } + + // create unfolded rec points and fill them with new energy lists + // First calculate energy deposited in each sell in accordance with + // fit (without fluctuations): efit[] + // and later correct this number in acordance with actual energy + // deposition + + Int_t nDigits = iniTower->GetMultiplicity() ; + Float_t * efit = new Float_t[nDigits] ; + Double_t xDigit=0.,yDigit=0.,zDigit=0. ; + Float_t xpar=0.,zpar=0.,epar=0. ; + + AliEMCALDigit * digit = 0 ; + Int_t * Digits = iniTower->GetDigitsList() ; + + Int_t iparam ; + Int_t iDigit ; + for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ + digit = dynamic_cast( fDigitsArr->At(Digits[iDigit] ) ) ; + fGeom->RelPosCellInSModule(digit->GetId(), yDigit, xDigit, zDigit); + efit[iDigit] = 0; + + iparam = 0 ; + while(iparam < nPar ){ + xpar = fitparameters[iparam] ; + zpar = fitparameters[iparam+1] ; + epar = fitparameters[iparam+2] ; + iparam += 3 ; + efit[iDigit] += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ; + } + } + + + // Now create new RecPoints and fill energy lists with efit corrected to fluctuations + // so that energy deposited in each cell is distributed between new clusters proportionally + // to its contribution to efit + + Float_t * Energies = iniTower->GetEnergiesList() ; + Float_t ratio ; + + iparam = 0 ; + while(iparam < nPar ){ + xpar = fitparameters[iparam] ; + zpar = fitparameters[iparam+1] ; + epar = fitparameters[iparam+2] ; + iparam += 3 ; + + AliEMCALRecPoint * RecPoint = 0 ; + + if(fNumberOfECAClusters >= fRecPoints->GetSize()) + fRecPoints->Expand(2*fNumberOfECAClusters) ; + + (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ; + RecPoint = dynamic_cast( fRecPoints->At(fNumberOfECAClusters) ) ; + fNumberOfECAClusters++ ; + RecPoint->SetNExMax((Int_t)nPar/3) ; + + Float_t eDigit ; + for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ + digit = dynamic_cast( fDigitsArr->At( Digits[iDigit] ) ) ; + fGeom->RelPosCellInSModule(digit->GetId(), yDigit, xDigit, zDigit); + + ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ; + eDigit = Energies[iDigit] * ratio ; + RecPoint->AddDigit( *digit, eDigit ) ; + } + } + + delete[] fitparameters ; + delete[] efit ; } //_____________________________________________________________________________ -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 - - ::Fatal("UnfoldingChiSquare","Unfolding not implemented") ; + + TList * toMinuit = dynamic_cast( gMinuit->GetObjectFit() ) ; + + AliEMCALRecPoint * RecPoint = dynamic_cast( toMinuit->At(0) ) ; + TClonesArray * digits = dynamic_cast( toMinuit->At(1) ) ; + // A bit buggy way to get an access to the geometry + // To be revised! + AliEMCALGeometry *geom = dynamic_cast(toMinuit->At(2)); + + Int_t * Digits = RecPoint->GetDigitsList() ; + + Int_t nOdigits = RecPoint->GetDigitsMultiplicity() ; + + Float_t * Energies = RecPoint->GetEnergiesList() ; + + fret = 0. ; + Int_t iparam ; + + if(iflag == 2) + for(iparam = 0 ; iparam < nPar ; iparam++) + Grad[iparam] = 0 ; // Will evaluate gradient + + Double_t efit ; + + AliEMCALDigit * digit ; + Int_t iDigit ; + + for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) { + + digit = dynamic_cast( digits->At( Digits[iDigit] ) ); + + Double_t xDigit=0 ; + Double_t zDigit=0 ; + Double_t yDigit=0 ;//not used yet, assumed to be 0 + + geom->RelPosCellInSModule(digit->GetId(), yDigit, xDigit, zDigit); + + if(iflag == 2){ // calculate gradient + Int_t iParam = 0 ; + efit = 0 ; + while(iParam < nPar ){ + Double_t dx = (xDigit - x[iParam]) ; + iParam++ ; + Double_t dz = (zDigit - x[iParam]) ; + iParam++ ; + efit += x[iParam] * ShowerShape(dx,dz) ; + iParam++ ; + } + Double_t sum = 2. * (efit - Energies[iDigit]) / Energies[iDigit] ; // Here we assume, that sigma = sqrt(E) + iParam = 0 ; + while(iParam < nPar ){ + Double_t xpar = x[iParam] ; + Double_t zpar = x[iParam+1] ; + Double_t epar = x[iParam+2] ; + Double_t dr = TMath::Sqrt( (xDigit - xpar) * (xDigit - xpar) + (zDigit - zpar) * (zDigit - zpar) ); + Double_t shape = sum * ShowerShape(xDigit - xpar,zDigit - zpar) ; + Double_t r133 = TMath::Power(dr, 1.33); + Double_t r669 = TMath::Power(dr,6.69); + Double_t deriv =-1.33 * TMath::Power(dr,0.33)*dr * ( 1.57 / ( (1.57 + 0.0860 * r133) * (1.57 + 0.0860 * r133) ) + - 0.55 / (1 + 0.000563 * r669) / ( (1 + 0.000563 * r669) * (1 + 0.000563 * r669) ) ) ; + + Grad[iParam] += epar * shape * deriv * (xpar - xDigit) ; // Derivative over x + iParam++ ; + Grad[iParam] += epar * shape * deriv * (zpar - zDigit) ; // Derivative over z + iParam++ ; + Grad[iParam] += shape ; // Derivative over energy + iParam++ ; + } + } + efit = 0; + iparam = 0 ; + + + while(iparam < nPar ){ + Double_t xpar = x[iparam] ; + Double_t zpar = x[iparam+1] ; + Double_t epar = x[iparam+2] ; + iparam += 3 ; + efit += epar * ShowerShape(xDigit - xpar,zDigit - zpar) ; + } + + fret += (efit-Energies[iDigit])*(efit-Energies[iDigit])/Energies[iDigit] ; + // Here we assume, that sigma = sqrt(E) + } } //____________________________________________________________________________ void AliEMCALClusterizerv1::Print(Option_t * /*option*/)const @@ -788,18 +776,12 @@ void AliEMCALClusterizerv1::Print(Option_t * /*option*/)const // Print parameters - TString taskName(GetName()) ; - taskName.ReplaceAll(Version(), "") ; + TString taskName(Version()) ; printf("--------------- "); printf(taskName.Data()) ; printf(" "); - printf(GetTitle()) ; - printf("-----------\n"); - printf("Clusterizing digits from the file: "); - printf(taskName.Data()); - printf("\n Branch: "); - printf(GetName()); + printf("Clusterizing digits: "); printf("\n ECA Local Maximum cut = %f", fECALocMaxCut); printf("\n ECA Logarithmic weight = %f", fECAW0); if(fToUnfold) @@ -817,19 +799,13 @@ void AliEMCALClusterizerv1::Print(Option_t * /*option*/)const void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option) { // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer - AliEMCALLoader *emcalLoader = dynamic_cast(AliRunLoader::GetRunLoader()->GetDetectorLoader("EMCAL")); - TObjArray * aECARecPoints = emcalLoader->RecPoints() ; - if(strstr(option,"deb")) { printf("PrintRecPoints: Clusterization result:") ; - printf("event # %d\n", emcalLoader->GetRunLoader()->GetEventNumber() ) ; printf(" Found %d ECA Rec Points\n ", - aECARecPoints->GetEntriesFast()) ; + fRecPoints->GetEntriesFast()) ; } - fRecPointsInRun += aECARecPoints->GetEntriesFast() ; - if(strstr(option,"all")) { if(strstr(option,"deb")) { printf("\n-----------------------------------------------------------------------\n") ; @@ -837,15 +813,9 @@ void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option) 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)) ; + for (index = 0 ; index < fRecPoints->GetEntries() ; index++) { + AliEMCALRecPoint * rp = dynamic_cast(fRecPoints->At(index)) ; TVector3 globalpos; //rp->GetGlobalPosition(globalpos); TVector3 localpos; @@ -860,19 +830,6 @@ void AliEMCALClusterizerv1::PrintRecPoints(Option_t * option) 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()); - ///////////// if(strstr(option,"deb")){ for (Int_t iprimary=0; iprimaryFill(maxE); - fMaxL1->Fill(maxL1); - fMaxL2->Fill(maxL2); - fMaxDis->Fill(maxDis); - if(strstr(option,"deb")) printf("\n-----------------------------------------------------------------------\n"); } } -TList* AliEMCALClusterizerv1::BookHists() -{ - //set up histograms for monitoring clusterizer performance - - gROOT->cd(); - - fPointE = new TH1F("00_pointE","point energy", 2000, 0.0, 150.); - fPointL1 = new TH1F("01_pointL1","point L1", 1000, 0.0, 3.); - fPointL2 = new TH1F("02_pointL2","point L2", 1000, 0.0, 3.); - fPointDis = new TH1F("03_pointDisp","point dispersion", 1000, 0.0, 10.); - fPointMult = new TH1F("04_pointMult","#cell in point(cluster)", 101, -0.5, 100.5); - fDigitAmp = new TH1F("05_digitAmp","Digit Amplitude", 2000, 0.0, 5000.); - fMaxE = new TH1F("06_maxE","Max point energy", 2000, 0.0, 150.); - fMaxL1 = new TH1F("07_maxL1","Largest (first) of eigenvalue of covariance matrix", 1000, 0.0, 3.); - fMaxL2 = new TH1F("08_maxL2","Smalest (second) of eigenvalue of covariace matrix", 1000, 0.0, 3.); - fMaxDis = new TH1F("09_maxDis","Point dispersion", 1000, 0.0, 10.); // 9 - // - new TH1F("10_adcOfDigits","adc of digits(threshold control)", 1001, -0.5, 1000.5); // 10 - new TH1F("11_energyOfDigits","energy of digits(threshold control)", 1000, 0.0, 1.); // 11 - new TH1F("12_numberOfPoints","number of points(clusters)", 101, -0.5, 100.5); // 12 - - return AliEMCALHistoUtilities::MoveHistsToList("EmcalClusterizerv1ControlHists", kFALSE); -} - -void AliEMCALClusterizerv1::SaveHists(const char *fn) -{ - AliEMCALHistoUtilities::SaveListOfHists(fHists, fn, kTRUE); -} +//___________________________________________________________________ void AliEMCALClusterizerv1::PrintRecoInfo() { printf(" AliEMCALClusterizerv1::PrintRecoInfo() : version %s \n", Version() ); - TH1F *h = (TH1F*)fHists->At(12); - if(h) { - printf(" ## Multiplicity of RecPoints ## \n"); - for(int i=1; i<=h->GetNbinsX(); i++) { - int nbin = int((*h)[i]); - int mult = int(h->GetBinCenter(i)); - if(nbin > 0) printf(" %i : %5.5i %6.3f %% \n", mult, nbin, 100.*nbin/h->GetEntries()); - } - } -} -void AliEMCALClusterizerv1::DrawLambdasHists() -{ - if(fMaxL1) { - fMaxL1->Draw(); - if(fMaxL2) fMaxL2->Draw("same"); - if(fMaxDis) { - fMaxDis->Draw("same"); - } - } -} - -void AliEMCALClusterizerv1::Browse(TBrowser* b) -{ - if(fHists) b->Add(fHists); - if(fGeom) b->Add(fGeom); - TTask::Browse(b); }