X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliEMCALClusterizerNxN.cxx;h=5d0be0dca87f73f83db68e510af8789066b1a426;hb=53918c82f01363eb4ea57b702b513777b9234cad;hp=75c5d92aba68e3ec098ad4087948bc09086604cc;hpb=65bec413f88781dcfc968082568d35c70258687f;p=u%2Fmrichter%2FAliRoot.git diff --git a/EMCAL/AliEMCALClusterizerNxN.cxx b/EMCAL/AliEMCALClusterizerNxN.cxx index 75c5d92aba6..5d0be0dca87 100644 --- a/EMCAL/AliEMCALClusterizerNxN.cxx +++ b/EMCAL/AliEMCALClusterizerNxN.cxx @@ -65,103 +65,91 @@ ClassImp(AliEMCALClusterizerNxN) -Bool_t AliEMCALClusterizerNxN::fgkIsInputCalibrated = kFALSE; - //____________________________________________________________________________ AliEMCALClusterizerNxN::AliEMCALClusterizerNxN() - : AliEMCALClusterizer() + : AliEMCALClusterizer(), fNRowDiff(1), fNColDiff(1), fEnergyGrad(0) { // ctor with the indication of the file where header Tree and digits Tree are stored } //____________________________________________________________________________ AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry) - : AliEMCALClusterizer(geometry) + : AliEMCALClusterizer(geometry), fNRowDiff(1), fNColDiff(1), fEnergyGrad(0) { // ctor with the indication of the file where header Tree and digits Tree are stored // use this contructor to avoid usage of Init() which uses runloader // change needed by HLT - MP - } //____________________________________________________________________________ AliEMCALClusterizerNxN::AliEMCALClusterizerNxN(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped) -: AliEMCALClusterizer(geometry, calib, caloped) +: AliEMCALClusterizer(geometry, calib, caloped), fNRowDiff(1), fNColDiff(1), fEnergyGrad(0) { - // ctor, geometry and calibration are initialized elsewhere. - + // ctor, geometry and calibration are initialized elsewhere. } - //____________________________________________________________________________ - AliEMCALClusterizerNxN::~AliEMCALClusterizerNxN() +AliEMCALClusterizerNxN::~AliEMCALClusterizerNxN() { // dtor } - //____________________________________________________________________________ void AliEMCALClusterizerNxN::Digits2Clusters(Option_t * option) { // Steering method to perform clusterization for the current event // in AliEMCALLoader - if(strstr(option,"tim")) + if (strstr(option,"tim")) gBenchmark->Start("EMCALClusterizer"); - if(strstr(option,"print")) - Print("") ; + if (strstr(option,"print")) + Print(""); //Get calibration parameters from file or digitizer default values. - GetCalibrationParameters() ; + GetCalibrationParameters(); //Get dead channel map from file or digitizer default values. - GetCaloCalibPedestal() ; + GetCaloCalibPedestal(); - fNumberOfECAClusters = 0; + MakeClusters(); //only the real clusters - MakeClusters() ; //only the real clusters - - if(fToUnfold){ + if (fToUnfold) { fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr); fClusterUnfolding->MakeUnfolding(); } //Evaluate position, dispersion and other RecPoint properties for EC section - Int_t index ; - for(index = 0; index < fRecPoints->GetEntries(); index++) - { + for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) { AliEMCALRecPoint * rp = dynamic_cast(fRecPoints->At(index)); - if(rp){ - rp->EvalAll(fECAW0,fDigitsArr) ; + if (rp) { + rp->EvalAll(fECAW0,fDigitsArr,fJustClusters); AliDebug(5, Form("MAX INDEX %d ", rp->GetMaximalEnergyIndex())); //For each rec.point set the distance to the nearest bad crystal - rp->EvalDistanceToBadChannels(fCaloPed); + if (fCaloPed) + rp->EvalDistanceToBadChannels(fCaloPed); } } - fRecPoints->Sort() ; + fRecPoints->Sort(); - for(index = 0; index < fRecPoints->GetEntries(); index++) - { - AliEMCALRecPoint * rp = dynamic_cast(fRecPoints->At(index)); - if(rp){ - rp->SetIndexInList(index) ; - rp->Print(); + for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) { + AliEMCALRecPoint *rp = dynamic_cast(fRecPoints->At(index)); + if (rp) { + rp->SetIndexInList(index); } else AliFatal("RecPoint NULL!!"); } - fTreeR->Fill(); + if (fTreeR) + fTreeR->Fill(); - if(strstr(option,"deb") || strstr(option,"all")) - PrintRecPoints(option) ; + if (strstr(option,"deb") || strstr(option,"all")) + PrintRecPoints(option); AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast())); - fRecPoints->Delete(); - - if(strstr(option,"tim")){ + if (strstr(option,"tim")) { gBenchmark->Stop("EMCALClusterizer"); printf("Exec took %f seconds for Clusterizing", gBenchmark->GetCpuTime("EMCALClusterizer")); @@ -171,81 +159,90 @@ void AliEMCALClusterizerNxN::Digits2Clusters(Option_t * option) //____________________________________________________________________________ Int_t AliEMCALClusterizerNxN::AreNeighbours(AliEMCALDigit * d1, AliEMCALDigit * d2, Bool_t & shared) const { - // Gives the neighbourness of two digits = 0 are not neighbour ; continue searching - // = 1 are neighbour - // = 2 is in different SM; continue searching - // In case it is in different SM, but same phi rack, check if neigbours at eta=0 - // neighbours are defined as digits having at least a common side - // 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 - - 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=0, coldiff=0; - - shared = kFALSE; - - fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1); - fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2); - fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1); - fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2); - - //If different SM, check if they are in the same phi, then consider cells close to eta=0 as neighbours; May 2010 - if(nSupMod1 != nSupMod2 ) - { - //Check if the 2 SM are in the same PHI position (0,1), (2,3), ... - Float_t smPhi1 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod1); - Float_t smPhi2 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod2); - - if(!TMath::AreEqualAbs(smPhi1, smPhi2, 1e-3)) return 2; //Not phi rack equal, not neighbours - - // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2 - // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0 - if(nSupMod1%2) ieta1+=AliEMCALGeoParams::fgkEMCALCols; - else ieta2+=AliEMCALGeoParams::fgkEMCALCols; - - shared = kTRUE; // maybe a shared cluster, we know this later, set it for the moment. - - }//Different SM, same phi - - rowdiff = TMath::Abs(iphi1 - iphi2); - coldiff = TMath::Abs(ieta1 - ieta2) ; + // Gives the neighbourness of two digits = 0 are not neighbour ; continue searching + // = 1 are neighbour + // = 2 is in different SM; continue searching + // In case it is in different SM, but same phi rack, check if neigbours at eta=0 + // neighbours are defined as digits having at least a common side + // 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 + + if (fEnergyGrad) { //false by default + if (d2->GetCalibAmp()>d1->GetCalibAmp()) + return 3; // energy of neighboring cell should be smaller in order to become a neighbor + } - // neighbours +-1 in col and row - if ( TMath::Abs(coldiff) < 2 && TMath::Abs(rowdiff) < 2) - { - - AliDebug(9, Form("AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n", - d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared)); - - return 1; - }//Neighbours - else - { - AliDebug(9, Form("NOT AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n", - d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared)); - shared = kFALSE; - return 2 ; - }//Not neighbours + Int_t nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0, iphi1=0, ieta1=0; + Int_t nSupMod2=0, nModule2=0, nIphi2=0, nIeta2=0, iphi2=0, ieta2=0; + Int_t rowdiff=0, coldiff=0; + + shared = kFALSE; + + fGeom->GetCellIndex(d1->GetId(), nSupMod1,nModule1,nIphi1,nIeta1); + fGeom->GetCellIndex(d2->GetId(), nSupMod2,nModule2,nIphi2,nIeta2); + fGeom->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, iphi1,ieta1); + fGeom->GetCellPhiEtaIndexInSModule(nSupMod2,nModule2,nIphi2,nIeta2, iphi2,ieta2); + + //If different SM, check if they are in the same phi, then consider cells close to eta=0 as neighbours; May 2010 + if (nSupMod1 != nSupMod2 ) + { + //Check if the 2 SM are in the same PHI position (0,1), (2,3), ... + Float_t smPhi1 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod1); + Float_t smPhi2 = fGeom->GetEMCGeometry()->GetPhiCenterOfSM(nSupMod2); + + if (!TMath::AreEqualAbs(smPhi1, smPhi2, 1e-3)) return 2; //Not phi rack equal, not neighbours + + // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2 + // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0 + if (nSupMod1%2) ieta1+=AliEMCALGeoParams::fgkEMCALCols; + else ieta2+=AliEMCALGeoParams::fgkEMCALCols; + + shared = kTRUE; // maybe a shared cluster, we know this later, set it for the moment. + + }//Different SM, same phi + + rowdiff = TMath::Abs(iphi1 - iphi2); + coldiff = TMath::Abs(ieta1 - ieta2); + + // neighbours +-1 in col and row + if ( TMath::Abs(coldiff) <= fNColDiff && TMath::Abs(rowdiff) <= fNRowDiff) + { + + AliDebug(9, Form("AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n", + d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared)); + + return 1; + }//Neighbours + else + { + AliDebug(9, Form("NOT AliEMCALClusterizerNxN::AreNeighbours(): id1=%d, (row %d, col %d) ; id2=%d, (row %d, col %d), shared %d \n", + d1->GetId(), iphi1,ieta1, d2->GetId(), iphi2,ieta2, shared)); + shared = kFALSE; + return 2; + }//Not neighbours } //____________________________________________________________________________ void AliEMCALClusterizerNxN::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 + // Make clusters - if (fGeom==0) AliFatal("Did not get geometry from EMCALLoader"); + if (fGeom==0) + AliFatal("Did not get geometry from EMCALLoader"); - fRecPoints->Clear(); + fNumberOfECAClusters = 0; + fRecPoints->Delete(); - // Set up TObjArray with pointers to digits to work on - //TObjArray *digitsC = new TObjArray(); + // Set up TObjArray with pointers to digits to work on, calibrate digits TObjArray digitsC; TIter nextdigit(fDigitsArr); AliEMCALDigit *digit = 0; - while ( (digit = dynamic_cast(nextdigit())) ) { + while ( (digit = static_cast(nextdigit())) ) { + Float_t dEnergyCalibrated = digit->GetAmplitude(); + Float_t time = digit->GetTime(); + Calibrate(dEnergyCalibrated,time ,digit->GetId()); + digit->SetCalibAmp(dEnergyCalibrated); + digit->SetTime(time); digitsC.AddLast(digit); } @@ -262,31 +259,27 @@ void AliEMCALClusterizerNxN::MakeClusters() Float_t dMaxEnergyDigit = -1; AliEMCALDigit *pMaxEnergyDigit = 0; nextdigitC.Reset(); - while ( (digit = dynamic_cast(nextdigitC())) ) + while ( (digit = static_cast(nextdigitC())) ) { // scan over the list of digitsC - Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId()); - //AliDebug(5, Form("-> Digit ENERGY: %1.5f", dEnergyCalibrated)); - - //if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > fECAClusteringThreshold ) - if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > 0.0) // no threshold! - { - if (dEnergyCalibrated > dMaxEnergyDigit) + Float_t dEnergyCalibrated = digit->GetCalibAmp(); + + if (fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > fMinECut) // no threshold by default! + { // needs to be set in OCDB! + if (dEnergyCalibrated > dMaxEnergyDigit) { dMaxEnergyDigit = dEnergyCalibrated; iMaxEnergyDigit = digit->GetId(); pMaxEnergyDigit = digit; } - } + } } - if (iMaxEnergyDigit < 0 || digitsC.GetEntries() <= 0) { bDone = kTRUE; continue; } - AliDebug (2, Form("Max digit found: %1.2f AbsId: %d", dMaxEnergyDigit, iMaxEnergyDigit)); - AliDebug(5, Form("Max Digit ENERGY: %1.5f", dMaxEnergyDigit)); + AliDebug (2, Form("Max digit found: %1.5f AbsId: %d", dMaxEnergyDigit, iMaxEnergyDigit)); // keep the candidate digits in a list TList clusterDigitList; @@ -295,51 +288,50 @@ void AliEMCALClusterizerNxN::MakeClusters() Double_t clusterCandidateEnergy = dMaxEnergyDigit; - // now loop over the resto of the digits and cluster into NxN cluster + // now loop over the rest of the digits and cluster into NxN cluster // we do not actually cluster yet: we keep them in the list clusterDigitList nextdigitC.Reset(); - while ( (digit = dynamic_cast(nextdigitC())) ) + while ( (digit = static_cast(nextdigitC())) ) { // scan over the list of digitsC if (digit == pMaxEnergyDigit) continue; - Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId()); + Float_t dEnergyCalibrated = digit->GetCalibAmp(); AliDebug(5, Form("-> Digit ENERGY: %1.5f", dEnergyCalibrated)); - //if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > fECAClusteringThreshold ) - if(fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > 0.0 ) - { - Float_t time = pMaxEnergyDigit->GetTime(); //Time or TimeR? - if(TMath::Abs(time - digit->GetTime()) > fTimeCut ) continue; //Time or TimeR? - Bool_t shared = kFALSE; //cluster shared by 2 SuperModules? - if (AreNeighbours(pMaxEnergyDigit, digit, shared) == 1) // call (digit,digitN) in THAT order !!!!! + if (fGeom->CheckAbsCellId(digit->GetId()) && dEnergyCalibrated > 0.0 ) + { + Float_t time = pMaxEnergyDigit->GetTime(); //Time or TimeR? + if (TMath::Abs(time - digit->GetTime()) > fTimeCut ) continue; //Time or TimeR? + Bool_t shared = kFALSE; //cluster shared by 2 SuperModules? + if (AreNeighbours(pMaxEnergyDigit, digit, shared) == 1) // call (digit,digitN) in THAT order !!!!! { - clusterDigitList.AddLast(digit) ; + clusterDigitList.AddLast(digit); clusterCandidateEnergy += dEnergyCalibrated; } - } + } }// loop over the next digits // start a cluster here only if a cluster energy is larger than clustering threshold - //if (clusterCandidateEnergy > 0.1) AliDebug(5, Form("Clusterization threshold is %f MeV", fECAClusteringThreshold)); if (clusterCandidateEnergy > fECAClusteringThreshold) { - if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(2*fNumberOfECAClusters+1) ; + if (fNumberOfECAClusters >= fRecPoints->GetSize()) + fRecPoints->Expand(2*fNumberOfECAClusters+1); - AliEMCALRecPoint *recPoint = new AliEMCALRecPoint("") ; - fRecPoints->AddAt(recPoint, fNumberOfECAClusters) ; - recPoint = dynamic_cast(fRecPoints->At(fNumberOfECAClusters)) ; - if(recPoint){ - fNumberOfECAClusters++ ; + (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("") ; + AliEMCALRecPoint *recPoint = dynamic_cast( fRecPoints->At(fNumberOfECAClusters) ) ; + // AliEMCALRecPoint *recPoint = new AliEMCALRecPoint(""); + // fRecPoints->AddAt(recPoint, fNumberOfECAClusters); + // recPoint = static_cast(fRecPoints->At(fNumberOfECAClusters)); + if (recPoint) { + fNumberOfECAClusters++; recPoint->SetClusterType(AliVCluster::kEMCALClusterv1); AliDebug(9, Form("Number of cells per cluster (max is 9!): %d", clusterDigitList.GetEntries())); for (Int_t idig = 0; idig < clusterDigitList.GetEntries(); idig++) { - digit = (AliEMCALDigit*)clusterDigitList.At(idig); - Float_t dEnergyCalibrated = Calibrate(digit->GetAmplitude(), digit->GetTime(),digit->GetId()); AliDebug(5, Form(" Adding digit %d", digit->GetId())); // note: this way the sharing info is lost! - recPoint->AddDigit(*digit, dEnergyCalibrated, kFALSE) ; //Time or TimeR? + recPoint->AddDigit(*digit, digit->GetCalibAmp(), kFALSE); //Time or TimeR? digitsC.Remove(digit); } }// recpoint @@ -355,16 +347,5 @@ void AliEMCALClusterizerNxN::MakeClusters() AliDebug (2, Form("Number of digits left: %d", digitsC.GetEntries())); } // while ! done - //delete digitsC ; //nope we use an object - AliDebug(1,Form("total no of clusters %d from %d digits",fNumberOfECAClusters,fDigitsArr->GetEntriesFast())); } - -//___________________________________________________________________ -void AliEMCALClusterizerNxN::SetInputCalibrated(Bool_t val) -{ - // - // input is calibrated - the case when we run already on ESD - // - AliEMCALClusterizerNxN::fgkIsInputCalibrated = val; -}