X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=PHOS%2FAliPHOSClusterizerv1.cxx;h=840ffd05454774e42831d4c4172e7c48c6b2a4ed;hb=04de96996745b8b2665d133da7836c302dbe028a;hp=6fb8aa30b62ea30c090a42559c50c630b5970aaa;hpb=6483babc2388acdec73dbabd46fa083362a36f58;p=u%2Fmrichter%2FAliRoot.git diff --git a/PHOS/AliPHOSClusterizerv1.cxx b/PHOS/AliPHOSClusterizerv1.cxx index 6fb8aa30b62..840ffd05454 100644 --- a/PHOS/AliPHOSClusterizerv1.cxx +++ b/PHOS/AliPHOSClusterizerv1.cxx @@ -180,7 +180,6 @@ #include "AliPHOSCpvRecPoint.h" #include "AliPHOSDigit.h" #include "AliPHOSDigitizer.h" -#include "AliPHOSCalibrationDB.h" #include "AliCDBManager.h" #include "AliCDBStorage.h" #include "AliCDBEntry.h" @@ -198,12 +197,17 @@ AliPHOSClusterizerv1::AliPHOSClusterizerv1() : fNumberOfEmcClusters(0), fNumberOfCpvClusters(0), fEmcClusteringThreshold(0), fCpvClusteringThreshold(0), fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0), - fW0CPV(0), fEmcTimeGate(0) + fW0CPV(0), + fTimeGateLowAmp(0.), fTimeGateLow(0.), fTimeGateHigh(0.), + fEcoreRadius(0) { // default ctor (to be used mainly by Streamer) - InitParameters() ; - fDefaultInit = kTRUE ; + fDefaultInit = kTRUE ; + + for(Int_t i=0; i<53760; i++){ + fDigitsUsed[i]=0 ; + } } //____________________________________________________________________________ @@ -214,11 +218,16 @@ AliPHOSClusterizerv1::AliPHOSClusterizerv1(AliPHOSGeometry *geom) : fNumberOfEmcClusters(0), fNumberOfCpvClusters(0), fEmcClusteringThreshold(0), fCpvClusteringThreshold(0), fEmcLocMaxCut(0), fW0(0), fCpvLocMaxCut(0), - fW0CPV(0), fEmcTimeGate(0) + fW0CPV(0), + fTimeGateLowAmp(0.), fTimeGateLow(0.), fTimeGateHigh(0.), + fEcoreRadius(0) { // ctor with the indication of the file where header Tree and digits Tree are stored - InitParameters() ; + for(Int_t i=0; i<53760; i++){ + fDigitsUsed[i]=0 ; + } + Init() ; fDefaultInit = kFALSE ; } @@ -277,6 +286,8 @@ Bool_t AliPHOSClusterizerv1::FindFit(AliPHOSEmcRecPoint * emcRP, AliPHOSDigit ** // Cluster will be fitted as a superposition of nPar/3 electromagnetic showers + if(!gMinuit) //it was deleted by someone else + gMinuit = new TMinuit(100) ; gMinuit->mncler(); // Reset Minuit's list of paramters gMinuit->SetPrintLevel(-1) ; // No Printout gMinuit->SetFCN(AliPHOSClusterizerv1::UnfoldingChiSquare) ; @@ -373,6 +384,8 @@ void AliPHOSClusterizerv1::Init() fgCalibData = new AliPHOSCalibData(-1); //use AliCDBManager's run number if (fgCalibData->GetCalibDataEmc() == 0) AliFatal("Calibration parameters for PHOS EMC not found. Stop reconstruction.\n"); + if (fgCalibData->GetCalibDataCpv() == 0) + AliFatal("Calibration parameters for PHOS CPV not found. Stop reconstruction.\n"); } @@ -383,24 +396,27 @@ void AliPHOSClusterizerv1::InitParameters() fNumberOfCpvClusters = 0 ; fNumberOfEmcClusters = 0 ; - const AliPHOSRecoParam* parEmc = AliPHOSReconstructor::GetRecoParamEmc(); - if(!parEmc) AliFatal("Reconstruction parameters for EMC not set!"); + const AliPHOSRecoParam* recoParam = AliPHOSReconstructor::GetRecoParam(); + if(!recoParam) AliFatal("Reconstruction parameters are not set!"); - const AliPHOSRecoParam* parCpv = AliPHOSReconstructor::GetRecoParamCpv(); - if(!parCpv) AliFatal("Reconstruction parameters for CPV not set!"); + recoParam->Print(); - fCpvClusteringThreshold = parCpv->GetClusteringThreshold(); - fEmcClusteringThreshold = parEmc->GetClusteringThreshold(); + fEmcClusteringThreshold = recoParam->GetEMCClusteringThreshold(); + fCpvClusteringThreshold = recoParam->GetCPVClusteringThreshold(); - fEmcLocMaxCut = parEmc->GetLocalMaxCut(); - fCpvLocMaxCut = parCpv->GetLocalMaxCut(); + fEmcLocMaxCut = recoParam->GetEMCLocalMaxCut(); + fCpvLocMaxCut = recoParam->GetCPVLocalMaxCut(); + + fW0 = recoParam->GetEMCLogWeight(); + fW0CPV = recoParam->GetCPVLogWeight(); - fW0 = parEmc->GetLogWeight(); - fW0CPV = parCpv->GetLogWeight(); + fTimeGateLowAmp = recoParam->GetTimeGateAmpThresh() ; + fTimeGateLow = recoParam->GetTimeGateLow() ; + fTimeGateHigh = recoParam->GetTimeGateHigh() ; - fEmcTimeGate = 1.e-6 ; + fEcoreRadius = recoParam->GetEMCEcoreRadius(); - fToUnfold = parEmc->ToUnfold() ; + fToUnfold = recoParam->EMCToUnfold() ; fWrite = kTRUE ; } @@ -428,7 +444,7 @@ Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)c if (( coldiff <= 1 ) && ( rowdiff <= 1 )){ //At least common vertex // if (( relid1[2]==relid2[2] && coldiff <= 1 ) || ( relid1[3]==relid2[3] && rowdiff <= 1 )){ //common side - if((relid1[1] != 0) || (TMath::Abs(d1->GetTime() - d2->GetTime() ) < fEmcTimeGate)) + if((relid1[1] != 0) || CheckTimeGate(d1->GetTime(),d1->GetEnergy(),d2->GetTime(),d2->GetEnergy())) return 1 ; } else { @@ -451,6 +467,18 @@ Int_t AliPHOSClusterizerv1::AreNeighbours(AliPHOSDigit * d1, AliPHOSDigit * d2)c return 0 ; } //____________________________________________________________________________ +Bool_t AliPHOSClusterizerv1::CheckTimeGate(Float_t t1, Float_t amp1, Float_t t2, Float_t amp2)const{ + //Check if two cells have reasonable time difference + //Note that at low amplitude time is defined up to 1 tick == 100 ns. + if(amp1GetEntriesFast(); - Float_t emcMinE= AliPHOSReconstructor::GetRecoParamEmc()->GetMinE(); //Minimal digit energy + Float_t emcMinE= AliPHOSReconstructor::GetRecoParam()->GetEMCMinE(); //Minimal digit energy + TVector3 fakeVtx(0.,0.,0.) ; for(index = 0; index < nEmc; index++){ AliPHOSEmcRecPoint * rp = - dynamic_cast( fEMCRecPoints->At(index) ); - rp->Purify(emcMinE) ; + static_cast( fEMCRecPoints->At(index) ); + rp->Purify(emcMinE,fDigitsArr) ; if(rp->GetMultiplicity()==0){ fEMCRecPoints->RemoveAt(index) ; delete rp ; @@ -500,8 +529,8 @@ void AliPHOSClusterizerv1::WriteRecPoints() } // No vertex is available now, calculate corrections in PID - rp->EvalAll(fW0,fDigitsArr) ; - TVector3 fakeVtx(0.,0.,0.) ; + rp->EvalAll(fDigitsArr) ; + rp->EvalCoreEnergy(fW0,fEcoreRadius,fDigitsArr) ; rp->EvalAll(fW0,fakeVtx,fDigitsArr) ; rp->EvalLocal2TrackingCSTransform(); } @@ -509,7 +538,7 @@ void AliPHOSClusterizerv1::WriteRecPoints() fEMCRecPoints->Sort() ; // fEMCRecPoints->Expand(fEMCRecPoints->GetEntriesFast()) ; for(index = 0; index < fEMCRecPoints->GetEntries(); index++){ - dynamic_cast( fEMCRecPoints->At(index) )->SetIndexInList(index) ; + static_cast( fEMCRecPoints->At(index) )->SetIndexInList(index) ; } //For each rec.point set the distance to the nearest bad crystal (BVP) @@ -517,14 +546,15 @@ void AliPHOSClusterizerv1::WriteRecPoints() //Now the same for CPV for(index = 0; index < fCPVRecPoints->GetEntries(); index++){ - AliPHOSCpvRecPoint * rp = dynamic_cast( fCPVRecPoints->At(index) ); - rp->EvalAll(fW0CPV,fDigitsArr) ; + AliPHOSCpvRecPoint * rp = static_cast( fCPVRecPoints->At(index) ); + rp->EvalAll(fDigitsArr) ; + rp->EvalAll(fW0CPV,fakeVtx,fDigitsArr) ; rp->EvalLocal2TrackingCSTransform(); } fCPVRecPoints->Sort() ; for(index = 0; index < fCPVRecPoints->GetEntries(); index++) - dynamic_cast( fCPVRecPoints->At(index) )->SetIndexInList(index) ; + static_cast( fCPVRecPoints->At(index) )->SetIndexInList(index) ; fCPVRecPoints->Expand(fCPVRecPoints->GetEntriesFast()) ; @@ -543,14 +573,16 @@ void AliPHOSClusterizerv1::MakeClusters() fNumberOfEmcClusters = 0 ; //Mark all digits as unused yet + const Int_t maxNDigits = 3584; // There is no clusters larger than PHOS module ;) Int_t nDigits=fDigitsArr->GetEntriesFast() ; + for(Int_t i=0; iGetEnergy() > fEmcClusteringThreshold ) || - ( IsInCpv(digit) && digit->GetEnergy() > fCpvClusteringThreshold ) ) { + if (( IsInEmc(digit) && Calibrate(digit->GetEnergy(),digit->GetId()) > fEmcClusteringThreshold ) || + ( IsInCpv(digit) && Calibrate(digit->GetEnergy(),digit->GetId()) > fCpvClusteringThreshold ) ) { Int_t iDigitInCluster = 0 ; - if ( IsInEmc(digit) ) { // start a new EMC RecPoint if(fNumberOfEmcClusters >= fEMCRecPoints->GetSize()) @@ -575,7 +606,7 @@ void AliPHOSClusterizerv1::MakeClusters() fEMCRecPoints->AddAt(new AliPHOSEmcRecPoint(""), fNumberOfEmcClusters) ; clu = static_cast( fEMCRecPoints->At(fNumberOfEmcClusters) ) ; fNumberOfEmcClusters++ ; - clu->AddDigit(*digit, digit->GetEnergy()) ; + clu->AddDigit(*digit, Calibrate(digit->GetEnergy(),digit->GetId()),CalibrateT(digit->GetTime(),digit->GetId(),digit->IsLG())) ; clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ; iDigitInCluster++ ; fDigitsUsed[i]=kTRUE ; @@ -587,12 +618,12 @@ void AliPHOSClusterizerv1::MakeClusters() fCPVRecPoints->AddAt(new AliPHOSCpvRecPoint(""), fNumberOfCpvClusters) ; clu = static_cast( fCPVRecPoints->At(fNumberOfCpvClusters) ) ; fNumberOfCpvClusters++ ; - clu->AddDigit(*digit, digit->GetEnergy()) ; + clu->AddDigit(*digit, Calibrate(digit->GetEnergy(),digit->GetId()),0.) ; // no timing information in CPV clusterdigitslist[iDigitInCluster] = digit->GetIndexInList() ; iDigitInCluster++ ; fDigitsUsed[i]=kTRUE ; } // else - + //Now scan remaining digits in list to find neigbours of our seed AliPHOSDigit * digitN ; @@ -601,6 +632,11 @@ void AliPHOSClusterizerv1::MakeClusters() digit = static_cast( fDigitsArr->At(clusterdigitslist[index]) ) ; index++ ; for(Int_t j=iFirst; j= maxNDigits) { + AliError(Form("The number of digits in cluster is more than %d, skip the rest of event", + maxNDigits)); + return; + } if(fDigitsUsed[j]) continue ; //look through remaining digits digitN = static_cast( fDigitsArr->At(j) ) ; @@ -612,22 +648,22 @@ void AliPHOSClusterizerv1::MakeClusters() case 0 : // not a neighbour break ; case 1 : // are neighbours - clu->AddDigit(*digitN, digitN->GetEnergy()); + clu->AddDigit(*digitN, Calibrate(digitN->GetEnergy(),digitN->GetId()),CalibrateT(digitN->GetTime(),digitN->GetId(),digit->IsLG())) ; clusterdigitslist[iDigitInCluster] = j ; iDigitInCluster++ ; fDigitsUsed[j]=kTRUE ; break ; case 2 : // too far from each other - goto endofloop; + goto endOfLoop; } // switch } - endofloop: ; //scanned all possible neighbours for this digit + endOfLoop: ; //scanned all possible neighbours for this digit } // loop over cluster } // energy theshold - } + } } @@ -646,7 +682,7 @@ void AliPHOSClusterizerv1::MakeUnfolding() Int_t index ; for(index = 0 ; index < numberofNotUnfolded ; index++){ - AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast( fEMCRecPoints->At(index) ) ; + AliPHOSEmcRecPoint * emcRecPoint = static_cast( fEMCRecPoints->At(index) ) ; if(emcRecPoint->GetPHOSMod()> nModulesToUnfold) break ; @@ -684,12 +720,12 @@ void AliPHOSClusterizerv1::MakeUnfolding() Int_t index ; for(index = 0 ; index < numberofCpvNotUnfolded ; index++){ - AliPHOSRecPoint * recPoint = dynamic_cast( fCPVRecPoints->At(index) ) ; + AliPHOSRecPoint * recPoint = static_cast( fCPVRecPoints->At(index) ) ; if(recPoint->GetPHOSMod()> nModulesToUnfold) break ; - AliPHOSEmcRecPoint * emcRecPoint = dynamic_cast(recPoint) ; + AliPHOSEmcRecPoint * emcRecPoint = static_cast(recPoint) ; Int_t nMultipl = emcRecPoint->GetMultiplicity() ; AliPHOSDigit ** maxAt = new AliPHOSDigit*[nMultipl] ; @@ -765,7 +801,7 @@ void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc, Int_t iparam ; Int_t iDigit ; for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ - digit = dynamic_cast( fDigitsArr->At(emcDigits[iDigit] ) ) ; + digit = static_cast( fDigitsArr->At(emcDigits[iDigit] ) ) ; fGeom->AbsToRelNumbering(digit->GetId(), relid) ; fGeom->RelPosInModule(relid, xDigit, zDigit) ; efit[iDigit] = 0; @@ -805,7 +841,7 @@ void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc, fEMCRecPoints->Expand(2*fNumberOfEmcClusters) ; (*fEMCRecPoints)[fNumberOfEmcClusters] = new AliPHOSEmcRecPoint("") ; - emcRP = dynamic_cast( fEMCRecPoints->At(fNumberOfEmcClusters) ) ; + emcRP = static_cast( fEMCRecPoints->At(fNumberOfEmcClusters) ) ; fNumberOfEmcClusters++ ; emcRP->SetNExMax((Int_t)nPar/3) ; } @@ -814,19 +850,19 @@ void AliPHOSClusterizerv1::UnfoldCluster(AliPHOSEmcRecPoint * iniEmc, fCPVRecPoints->Expand(2*fNumberOfCpvClusters) ; (*fCPVRecPoints)[fNumberOfCpvClusters] = new AliPHOSCpvRecPoint("") ; - emcRP = dynamic_cast( fCPVRecPoints->At(fNumberOfCpvClusters) ) ; + emcRP = static_cast( fCPVRecPoints->At(fNumberOfCpvClusters) ) ; fNumberOfCpvClusters++ ; } Float_t eDigit ; for(iDigit = 0 ; iDigit < nDigits ; iDigit ++){ - digit = dynamic_cast( fDigitsArr->At( emcDigits[iDigit] ) ) ; + digit = static_cast( fDigitsArr->At( emcDigits[iDigit] ) ) ; fGeom->AbsToRelNumbering(digit->GetId(), relid) ; fGeom->RelPosInModule(relid, xDigit, zDigit) ; // ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar,vIncid) / efit[iDigit] ; ratio = epar * ShowerShape(xDigit - xpar,zDigit - zpar) / efit[iDigit] ; eDigit = emcEnergies[iDigit] * ratio ; - emcRP->AddDigit( *digit, eDigit ) ; + emcRP->AddDigit( *digit, eDigit,CalibrateT(digit->GetTime(),digit->GetId(),digit->IsLG()) ) ; } } @@ -841,17 +877,17 @@ void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Dou // Calculates the Chi square for the cluster unfolding minimization // Number of parameters, Gradient, Chi squared, parameters, what to do - TList * toMinuit = dynamic_cast( gMinuit->GetObjectFit() ) ; + TList * toMinuit = static_cast( gMinuit->GetObjectFit() ) ; - AliPHOSEmcRecPoint * emcRP = dynamic_cast( toMinuit->At(0) ) ; - TClonesArray * digits = dynamic_cast( toMinuit->At(1) ) ; + AliPHOSEmcRecPoint * emcRP = static_cast( toMinuit->At(0) ) ; + TClonesArray * digits = static_cast( toMinuit->At(1) ) ; // A bit buggy way to get an access to the geometry // To be revised! - AliPHOSGeometry *geom = dynamic_cast(toMinuit->At(2)); + AliPHOSGeometry *geom = static_cast(toMinuit->At(2)); -// TVector3 * vtx = dynamic_cast(toMinuit->At(3)) ; //Vertex position +// TVector3 * vtx = static_cast(toMinuit->At(3)) ; //Vertex position - // AliPHOSEmcRecPoint * emcRP = dynamic_cast( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit + // AliPHOSEmcRecPoint * emcRP = static_cast( gMinuit->GetObjectFit() ) ; // EmcRecPoint to fit Int_t * emcDigits = emcRP->GetDigitsList() ; @@ -874,7 +910,7 @@ void AliPHOSClusterizerv1::UnfoldingChiSquare(Int_t & nPar, Double_t * Grad, Dou for( iDigit = 0 ; iDigit < nOdigits ; iDigit++) { - digit = dynamic_cast( digits->At( emcDigits[iDigit] ) ); + digit = static_cast( digits->At( emcDigits[iDigit] ) ); Int_t relid[4] ; Float_t xDigit ; @@ -970,7 +1006,7 @@ void AliPHOSClusterizerv1::Print(const Option_t *)const else message = " AliPHOSClusterizerv1 not initialized " ; - AliInfo(Form("%s, %s %s %s %s %s %s %s %s %s %s", message.Data(), + AliInfo(Form("%s, %s %s %s %s %f %f %f %f %f %f", message.Data(), taskName.Data(), GetTitle(), taskName.Data(), @@ -1038,9 +1074,9 @@ void AliPHOSClusterizerv1::SetDistancesToBadChannels() //Author: Boris Polichtchouk if(!fgCalibData->GetNumOfEmcBadChannels()) return; - AliInfo(Form("%d bad channel(s) found.\n",fgCalibData->GetNumOfEmcBadChannels())); Int_t badIds[8000]; + memset(badIds,0,8000*sizeof(Int_t)); fgCalibData->EmcBadChannelIds(badIds); AliPHOSEmcRecPoint* rp; @@ -1081,3 +1117,48 @@ void AliPHOSClusterizerv1::SetDistancesToBadChannels() } } +//================================================================================== +Float_t AliPHOSClusterizerv1::Calibrate(Float_t amp, Int_t absId) const{ + // Calibrate EMC digit, i.e. multiply its Amp by a factor read from CDB + + const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ; + + //Determine rel.position of the cell absolute ID + Int_t relId[4]; + geom->AbsToRelNumbering(absId,relId); + Int_t module=relId[0]; + Int_t row =relId[2]; + Int_t column=relId[3]; + if(relId[1]){ //CPV + Float_t calibration = fgCalibData->GetADCchannelCpv(module,column,row); + return amp*calibration ; + } + else{ //EMC + Float_t calibration = fgCalibData->GetADCchannelEmc(module,column,row); + return amp*calibration ; + } +} +//================================================================================== +Float_t AliPHOSClusterizerv1::CalibrateT(Float_t time, Int_t absId,Bool_t isLG)const{ + // Calibrate time in EMC digit + + const AliPHOSGeometry *geom = AliPHOSGeometry::GetInstance() ; + + //Determine rel.position of the cell absolute ID + Int_t relId[4]; + geom->AbsToRelNumbering(absId,relId); + Int_t module=relId[0]; + Int_t row =relId[2]; + Int_t column=relId[3]; + if(relId[1]){ //CPV + return 0. ; + } + else{ //EMC + if(isLG) + time += fgCalibData->GetLGTimeShiftEmc(module,column,row); + else + time += fgCalibData->GetTimeShiftEmc(module,column,row); + return time ; + } +} +