X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliEMCALRecPoint.cxx;h=a752edacda4ee416bbc5dd993f00dcdffd0b7985;hb=6c1053a8e7ff6e8fd16ddc33fb97fc4c6cf717d4;hp=96f08930a77a589350c9a2e3dd069578e38f4500;hpb=4800667c40c71eec8438d954e9b50ce50b05ebf5;p=u%2Fmrichter%2FAliRoot.git diff --git a/EMCAL/AliEMCALRecPoint.cxx b/EMCAL/AliEMCALRecPoint.cxx index 96f08930a77..a752edacda4 100644 --- a/EMCAL/AliEMCALRecPoint.cxx +++ b/EMCAL/AliEMCALRecPoint.cxx @@ -16,69 +16,124 @@ //_________________________________________________________________________ // Reconstructed Points for the EMCAL // A RecPoint is a cluster of digits +// +// //*-- Author: Yves Schutz (SUBATECH) //*-- Author: Dmitri Peressounko (RRC KI & SUBATECH) //*-- Author: Heather Gray (LBL) merged AliEMCALRecPoint and AliEMCALTowerRecPoint 02/04 // --- ROOT system --- -#include -#include -#include -#include -#include -#include +#include "TPad.h" +#include "TGraph.h" +#include "TPaveText.h" +#include "TClonesArray.h" +#include "TMath.h" +#include "TGeoMatrix.h" +#include "TGeoManager.h" +#include "TGeoPhysicalNode.h" // --- Standard library --- +#include // --- AliRoot header files --- -#include "AliGenerator.h" +//#include "AliGenerator.h" +class AliGenerator; +class AliEMCAL; +#include "AliLog.h" +#include "AliGeomManager.h" #include "AliEMCALGeometry.h" +#include "AliEMCALHit.h" #include "AliEMCALDigit.h" #include "AliEMCALRecPoint.h" +#include "AliCaloCalibPedestal.h" +#include "AliEMCALGeoParams.h" ClassImp(AliEMCALRecPoint) -AliEMCALGeometry * geom = 0; - //____________________________________________________________________________ AliEMCALRecPoint::AliEMCALRecPoint() - : AliRecPoint() + : AliCluster(), fGeomPtr(0), + fAmp(0), fIndexInList(-1), //to be set when the point is already stored + fLocPos(0,0,0), fLocPosM(0), + fMaxDigit(100), fMulDigit(0), fMaxTrack(200), + fMulTrack(0), fDigitsList(0), fTracksList(0), + fClusterType(-1), fCoreEnergy(0), fDispersion(0), + fEnergyList(0), fTimeList(0), fAbsIdList(0), + fTime(0.), fNExMax(0), fCoreRadius(10), //HG check this + fDETracksList(0), fMulParent(0), fMaxParent(0), + fParentsList(0), fDEParentsList(0), fSuperModuleNumber(0), + fDigitIndMax(-1), fDistToBadTower(-1) { // ctor - fMaxTrack = 0 ; - fMulDigit = 0 ; - fMaxParent = 0; - fMulParent = 0; - fAmp = 0. ; - fCoreEnergy = 0 ; - fEnergyList = 0 ; - fAbsIdList = 0; - fParentsList = 0; - fTime = 0. ; - // fLocPos.SetX(1.e+6) ; //Local position should be evaluated - fCoreRadius = 10; //HG Check this - geom = AliEMCALGeometry::GetInstance(); - geom->GetTransformationForSM(); // Global <-> Local + fGeomPtr = AliEMCALGeometry::GetInstance(); + + fLambda[0] = 0; + fLambda[1] = 0; + } //____________________________________________________________________________ -AliEMCALRecPoint::AliEMCALRecPoint(const char * opt) : AliRecPoint(opt) +AliEMCALRecPoint::AliEMCALRecPoint(const char *) + : AliCluster(), fGeomPtr(0), + fAmp(0), fIndexInList(-1), //to be set when the point is already stored + fLocPos(0,0,0), fLocPosM(new TMatrixF(3,3)), + fMaxDigit(100), fMulDigit(0), fMaxTrack(1000), fMulTrack(0), + fDigitsList(new Int_t[fMaxDigit]), fTracksList(new Int_t[fMaxTrack]), + fClusterType(-1), fCoreEnergy(0), fDispersion(0), + fEnergyList(new Float_t[fMaxDigit]), fTimeList(new Float_t[fMaxDigit]), + fAbsIdList(new Int_t[fMaxDigit]), fTime(-1.), fNExMax(0), fCoreRadius(10), + fDETracksList(new Float_t[fMaxTrack]), fMulParent(0), fMaxParent(1000), + fParentsList(new Int_t[fMaxParent]), fDEParentsList(new Float_t[fMaxParent]), + fSuperModuleNumber(0), fDigitIndMax(-1), fDistToBadTower(-1) { // ctor - fMaxTrack = 1000 ; - fMaxParent = 1000; - fMulDigit = 0 ; - fMulParent = 0; - fAmp = 0. ; - fCoreEnergy = 0 ; - fEnergyList = 0 ; - fAbsIdList = 0; - fParentsList = new Int_t[fMaxParent]; - fTime = -1. ; - //fLocPos.SetX(1.e+6) ; //Local position should be evaluated - fCoreRadius = 10; //HG Check this - geom = AliEMCALGeometry::GetInstance(); - geom->GetTransformationForSM(); // Global <-> Local + for (Int_t i = 0; i < fMaxTrack; i++) + fDETracksList[i] = 0; + for (Int_t i = 0; i < fMaxParent; i++) { + fParentsList[i] = -1; + fDEParentsList[i] = 0; + } + + fGeomPtr = AliEMCALGeometry::GetInstance(); + fLambda[0] = 0; + fLambda[1] = 0; +} + +//____________________________________________________________________________ +AliEMCALRecPoint::AliEMCALRecPoint(const AliEMCALRecPoint & rp) + : AliCluster(rp), fGeomPtr(rp.fGeomPtr), + fAmp(rp.fAmp), fIndexInList(rp.fIndexInList), + fLocPos(rp.fLocPos), fLocPosM(rp.fLocPosM), + fMaxDigit(rp.fMaxDigit), fMulDigit(rp.fMulDigit), + fMaxTrack(rp.fMaxTrack), fMulTrack(rp.fMaxTrack), + fDigitsList(new Int_t[rp.fMaxDigit]), fTracksList(new Int_t[rp.fMaxTrack]), + fClusterType(rp.fClusterType), fCoreEnergy(rp.fCoreEnergy), + fDispersion(rp.fDispersion), + fEnergyList(new Float_t[rp.fMaxDigit]), fTimeList(new Float_t[rp.fMaxDigit]), + fAbsIdList(new Int_t[rp.fMaxDigit]), fTime(rp.fTime), fNExMax(rp.fNExMax),fCoreRadius(rp.fCoreRadius), + fDETracksList(new Float_t[rp.fMaxTrack]), fMulParent(rp.fMulParent), + fMaxParent(rp.fMaxParent), fParentsList(new Int_t[rp.fMaxParent]), + fDEParentsList(new Float_t[rp.fMaxParent]), + fSuperModuleNumber(rp.fSuperModuleNumber), fDigitIndMax(rp.fDigitIndMax), + fDistToBadTower(rp.fDistToBadTower) +{ + //copy ctor + fLambda[0] = rp.fLambda[0]; + fLambda[1] = rp.fLambda[1]; + + for(Int_t i = 0; i < rp.fMulDigit; i++) { + fEnergyList[i] = rp.fEnergyList[i]; + fTimeList[i] = rp.fTimeList[i]; + fAbsIdList[i] = rp.fAbsIdList[i]; + } + + for(Int_t i = 0; i < rp.fMulTrack; i++) fDETracksList[i] = rp.fDETracksList[i]; + + for(Int_t i = 0; i < rp.fMulParent; i++) { + fParentsList[i] = rp.fParentsList[i]; + fDEParentsList[i] = rp.fDEParentsList[i]; + } + } //____________________________________________________________________________ AliEMCALRecPoint::~AliEMCALRecPoint() @@ -86,10 +141,68 @@ AliEMCALRecPoint::~AliEMCALRecPoint() // dtor if ( fEnergyList ) delete[] fEnergyList ; + if ( fTimeList ) + delete[] fTimeList ; if ( fAbsIdList ) delete[] fAbsIdList ; + if ( fDETracksList) + delete[] fDETracksList; if ( fParentsList) delete[] fParentsList; + if ( fDEParentsList) + delete[] fDEParentsList; + + delete fLocPosM ; + delete [] fDigitsList ; + delete [] fTracksList ; +} + +//____________________________________________________________________________ +AliEMCALRecPoint& AliEMCALRecPoint::operator= (const AliEMCALRecPoint &rp) +{ + // assignment operator + + if(&rp == this) return *this; + + fGeomPtr = rp.fGeomPtr; + fAmp = rp.fAmp; + fIndexInList = rp.fIndexInList; + fLocPos = rp.fLocPos; + fLocPosM = rp.fLocPosM; + fMaxDigit = rp.fMaxDigit; + fMulDigit = rp.fMulDigit; + fMaxTrack = rp.fMaxTrack; + fMulTrack = rp.fMaxTrack; + for(Int_t i = 0; iGetSuperModuleNumber(digit.GetId()); } if ( fMulDigit >= fMaxDigit ) { // increase the size of the lists fMaxDigit*=2 ; Int_t * tempo = new Int_t[fMaxDigit]; Float_t * tempoE = new Float_t[fMaxDigit]; + Float_t * tempoT = new Float_t[fMaxDigit]; Int_t * tempoId = new Int_t[fMaxDigit]; Int_t index ; for ( index = 0 ; index < fMulDigit ; index++ ){ tempo[index] = fDigitsList[index] ; tempoE[index] = fEnergyList[index] ; + tempoT[index] = fTimeList[index] ; tempoId[index] = fAbsIdList[index] ; } - delete [] fDigitsList ; - fDigitsList = new Int_t[fMaxDigit]; - + delete [] fDigitsList ; delete [] fEnergyList ; - fEnergyList = new Float_t[fMaxDigit]; - + delete [] fTimeList ; delete [] fAbsIdList ; - fAbsIdList = new Int_t[fMaxDigit]; - for ( index = 0 ; index < fMulDigit ; index++ ){ - fDigitsList[index] = tempo[index] ; - fEnergyList[index] = tempoE[index] ; - fAbsIdList[index] = tempoId[index] ; - } - - delete [] tempo ; - delete [] tempoE ; - delete [] tempoId ; + fDigitsList = tempo; + fEnergyList = tempoE; + fTimeList = tempoT; + fAbsIdList = tempoId; } // if fDigitsList[fMulDigit] = digit.GetIndexInList() ; fEnergyList[fMulDigit] = Energy ; + fTimeList[fMulDigit] = digit.GetTimeR() ; fAbsIdList[fMulDigit] = digit.GetId(); fMulDigit++ ; fAmp += Energy ; + //JLK 10-Oct-2007 this hasn't been filled before because it was in + //the wrong place in previous versions. + //Now we evaluate it only if the supermodulenumber for this recpoint + //has not yet been set (or is the 0th one) + if(fSuperModuleNumber == 0) + fSuperModuleNumber = fGeomPtr->GetSuperModuleNumber(digit.GetId()); + } //____________________________________________________________________________ Bool_t AliEMCALRecPoint::AreNeighbours(AliEMCALDigit * digit1, AliEMCALDigit * digit2 ) const @@ -152,20 +267,18 @@ Bool_t AliEMCALRecPoint::AreNeighbours(AliEMCALDigit * digit1, AliEMCALDigit * d // A neighbour is defined as being two digits which share a corner static Bool_t areNeighbours = kFALSE ; - static Int_t nSupMod=0, nTower=0, nIphi=0, nIeta=0; - static int nSupMod1=0, nTower1=0, nIphi1=0, nIeta1=0; + static Int_t nSupMod=0, nModule=0, nIphi=0, nIeta=0; + static int nSupMod1=0, nModule1=0, nIphi1=0, nIeta1=0; static Int_t relid1[2] , relid2[2] ; // ieta, iphi static Int_t rowdiff=0, coldiff=0; areNeighbours = kFALSE ; - // AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance(); - - geom->GetCellIndex(digit1->GetId(), nSupMod,nTower,nIphi,nIeta); - geom->GetCellPhiEtaIndexInSModule(nSupMod,nTower,nIphi,nIeta, relid1[0],relid1[1]); + fGeomPtr->GetCellIndex(digit1->GetId(), nSupMod,nModule,nIphi,nIeta); + fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, relid1[0],relid1[1]); - geom->GetCellIndex(digit2->GetId(), nSupMod1,nTower1,nIphi1,nIeta1); - geom->GetCellPhiEtaIndexInSModule(nSupMod1,nTower1,nIphi1,nIeta1, relid2[0],relid2[1]); + fGeomPtr->GetCellIndex(digit2->GetId(), nSupMod1,nModule1,nIphi1,nIeta1); + fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod1,nModule1,nIphi1,nIeta1, relid2[0],relid2[1]); rowdiff = TMath::Abs( relid1[0] - relid2[0] ) ; coldiff = TMath::Abs( relid1[1] - relid2[1] ) ; @@ -193,7 +306,7 @@ Int_t AliEMCALRecPoint::Compare(const TObject * obj) const TVector3 locpos2; clu->GetLocalPosition(locpos2); - Int_t rowdif = (Int_t)TMath::Ceil(locpos1.X()/delta)-(Int_t)TMath::Ceil(locpos2.X()/delta) ; + Int_t rowdif = (Int_t)(TMath::Ceil(locpos1.X()/delta)-TMath::Ceil(locpos2.X()/delta)) ; if (rowdif> 0) rv = 1 ; else if(rowdif < 0) @@ -246,7 +359,7 @@ void AliEMCALRecPoint::ExecuteEvent(Int_t /*event*/, Int_t, Int_t) // static Int_t pxold, pyold; - /* static TGraph * digitgraph = 0 ; + /* static TGraph * digitgraph = 0 ; static TPaveText* clustertext = 0 ; if (!gPad->IsEditable()) return; @@ -256,7 +369,6 @@ void AliEMCALRecPoint::ExecuteEvent(Int_t /*event*/, Int_t, Int_t) case kButton1Down:{ AliEMCALDigit * digit ; - AliEMCALGeometry * emcalgeom = (AliEMCALGetter::Instance())->EMCALGeometry() ; Int_t iDigit; Int_t relid[2] ; @@ -268,8 +380,8 @@ void AliEMCALRecPoint::ExecuteEvent(Int_t /*event*/, Int_t, Int_t) for(iDigit = 0; iDigit < kMulDigit; iDigit++) { Fatal("AliEMCALRecPoint::ExecuteEvent", " -> Something wrong with the code"); digit = 0 ; //dynamic_cast((fDigitsList)[iDigit]); - emcalgeom->AbsToRelNumbering(digit->GetId(), relid) ; - emcalgeom->PosInAlice(relid, xi[iDigit], zi[iDigit]) ; + fGeomPtr->AbsToRelNumbering(digit->GetId(), relid) ; + fGeomPtr->PosInAlice(relid, xi[iDigit], zi[iDigit]) ; } if (!digitgraph) { @@ -318,51 +430,70 @@ break; void AliEMCALRecPoint::EvalAll(Float_t logWeight,TClonesArray * digits) { // Evaluates all shower parameters - EvalLocalPosition(logWeight, digits) ; - // printf("eval position done\n"); EvalElipsAxis(logWeight, digits) ; - // printf("eval axis done\n"); EvalDispersion(logWeight, digits) ; - // printf("eval dispersion done\n"); - // EvalCoreEnergy(logWeight, digits); - // printf("eval energy done\n"); + //EvalCoreEnergy(logWeight, digits); EvalTime(digits) ; - // printf("eval time done\n"); - EvalPrimaries(digits) ; - // printf("eval pri done\n"); EvalParents(digits); - // printf("eval parent done\n"); + + //Called last because it sets the global position of the cluster? + EvalLocal2TrackingCSTransform(); + } //____________________________________________________________________________ void AliEMCALRecPoint::EvalDispersion(Float_t logWeight, TClonesArray * digits) { // Calculates the dispersion of the shower at the origin of the RecPoint + // in cell units - Nov 16,2006 - Double_t d = 0., wtot = 0., w = 0., xyzi[3], diff=0.; - Int_t iDigit=0, nstat=0, i=0; + Double_t d = 0., wtot = 0., w = 0.; + Int_t iDigit=0, nstat=0; AliEMCALDigit * digit ; - // AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance(); - - // Calculates the centre of gravity in the local EMCAL-module coordinates - if (!fLocPos.Mag()) - EvalLocalPosition(logWeight, digits) ; - // Calculates the dispersion in coordinates + // Calculates the dispersion in cell units + Double_t etai, phii, etaMean=0.0, phiMean=0.0; + int nSupMod=0, nModule=0, nIphi=0, nIeta=0; + int iphi=0, ieta=0; + // Calculate mean values for(iDigit=0; iDigit < fMulDigit; iDigit++) { digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit]) ; - geom->RelPosCellInSModule(digit->GetId(), xyzi[0], xyzi[1], xyzi[2]); - w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; + if (fAmp>0 && fEnergyList[iDigit]>0) { + fGeomPtr->GetCellIndex(digit->GetId(), nSupMod,nModule,nIphi,nIeta); + fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta); + etai=(Double_t)ieta; + phii=(Double_t)iphi; + w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; + + if(w>0.0) { + phiMean += phii*w; + etaMean += etai*w; + wtot += w; + } + } + } + if (wtot>0) { + phiMean /= wtot ; + etaMean /= wtot ; + } else AliError(Form("Wrong weight %f\n", wtot)); - if(w>0.0) { - wtot += w; - nstat++; - for(i=0; i<3; i++ ) { - diff = xyzi[i] - double(fLocPos[i]); - d += w * diff*diff; + // Calculate dispersion + for(iDigit=0; iDigit < fMulDigit; iDigit++) { + digit = (AliEMCALDigit *) digits->At(fDigitsList[iDigit]) ; + + if (fAmp>0 && fEnergyList[iDigit]>0) { + fGeomPtr->GetCellIndex(digit->GetId(), nSupMod,nModule,nIphi,nIeta); + fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta); + etai=(Double_t)ieta; + phii=(Double_t)iphi; + w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; + + if(w>0.0) { + nstat++; + d += w*((etai-etaMean)*(etai-etaMean)+(phii-phiMean)*(phii-phiMean)); } } } @@ -373,25 +504,77 @@ void AliEMCALRecPoint::EvalDispersion(Float_t logWeight, TClonesArray * digits) fDispersion = TMath::Sqrt(d) ; } +//____________________________________________________________________________ +void AliEMCALRecPoint::EvalDistanceToBadChannels(AliCaloCalibPedestal* caloped) +{ + //For each EMC rec. point set the distance to the nearest bad channel. + //AliInfo(Form("%d bad channel(s) found.\n", caloped->GetDeadTowerCount())); + //Needs to be carefully checked!!! Gustavo 10-11-2009 + + if(!caloped->GetDeadTowerCount()) return; + + //Number of supermodule this cluster belongs. + Int_t iSM = GetSuperModuleNumber(); + + //Get channels map of the supermodule where the cluster is + TH2D* hMap = caloped->GetDeadMap(iSM); + + TVector3 dR; + TVector3 cellpos; + + Float_t minDist = 100000; + Float_t dist = 0; + Int_t absId =-1; + + //Loop on tower status map + for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){ + for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){ + //Check if tower is bad. + if(hMap->GetBinContent(icol,irow)==AliCaloCalibPedestal::kAlive) continue; + //printf("AliEMCALRecPoint::EvalDistanceToBadChannels() - Bad channel in SM %d, col %d, row %d\n",iSM,icol, irow); + //Tower is bad, get the absId of the index. + absId = fGeomPtr->GetAbsCellIdFromCellIndexes(iSM, irow, icol); + //Get the position of this tower. + fGeomPtr->RelPosCellInSModule(absId,cellpos); + //Calculate distance between this tower and cluster, set if is smaller than previous. + dR = cellpos-fLocPos; + dist = dR.Mag(); + if(dist(digits->At(fDigitsList[iDigit])) ; + if(iDigit==0) { + idMax = digit->GetId(); // is it correct + dist = TmaxInCm(Double_t(fAmp)); + } + fGeomPtr->RelPosCellInSModule(digit->GetId(), idMax, dist, xyzi[0], xyzi[1], xyzi[2]); + //printf(" Id %i : dist %f Local x,y,z %f %f %f \n", digit->GetId(), dist, xyzi[0], xyzi[1], xyzi[2]); - geom->RelPosCellInSModule(digit->GetId(), xyzi[0], xyzi[1], xyzi[2]); - // printf(" Id %i : Local x,y,z %f %f %f \n", digit->GetId(), xyzi[0], xyzi[1], xyzi[2]); + //fGeomPtr->RelPosCellInSModule(digit->GetId(), xyzi[0], xyzi[1], xyzi[2]); + //printf(" Id %i : dist %f Local x,y,z %f %f %f \n", digit->GetId(), 0.0, xyzi[0], xyzi[1], xyzi[2]); + // if(fAmp>102.) assert(0); - if(logWeight > 0.0) w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp )); - else w = fEnergyList[iDigit]; // just energy + if(logWeight > 0.0) w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp )); + else w = fEnergyList[iDigit]; // just energy if(w>0.0) { wtot += w ; @@ -430,21 +613,181 @@ void AliEMCALRecPoint::EvalLocalPosition(Float_t logWeight, TClonesArray * digit fLocPosM = 0 ; // covariance matrix } -//void AliEMCALRecPoint::EvalLocalPositionSimple() -//{ // Weight is proportional of cell energy -//} +//____________________________________________________________________________ +void AliEMCALRecPoint::EvalLocalPositionFit(Double_t deff, Double_t logWeight, +Double_t phiSlope, TClonesArray * digits) +{ + // Aug 14-16, 2007 - for fit + // Aug 31 - should be static ?? + static Double_t dist, ycorr; + static AliEMCALDigit *digit; + + Int_t i=0, nstat=0, idMax=-1; + Double_t clXYZ[3]={0.,0.,0.}, clRmsXYZ[3]={0.,0.,0.}, xyzi[3], wtot=0., w=0.; + + for(Int_t iDigit=0; iDigitGetEntries(); iDigit++) { + digit = dynamic_cast(digits->At(fDigitsList[iDigit])) ; + if(iDigit==0) { + idMax = digit->GetId(); // is it correct + dist = TmaxInCm(Double_t(fAmp)); + } + + dist = deff; + fGeomPtr->RelPosCellInSModule(digit->GetId(), idMax, dist, xyzi[0], xyzi[1], xyzi[2]); + + if(logWeight > 0.0) w = TMath::Max( 0., logWeight + TMath::Log( fEnergyList[iDigit] / fAmp )); + else w = fEnergyList[iDigit]; // just energy + + if(w>0.0) { + wtot += w ; + nstat++; + for(i=0; i<3; i++ ) { + clXYZ[i] += (w*xyzi[i]); + clRmsXYZ[i] += (w*xyzi[i]*xyzi[i]); + } + } + } + // cout << " wtot " << wtot << endl; + if ( wtot > 0 ) { + // xRMS = TMath::Sqrt(x2m - xMean*xMean); + for(i=0; i<3; i++ ) { + clXYZ[i] /= wtot; + if(nstat>1) { + clRmsXYZ[i] /= (wtot*wtot); + clRmsXYZ[i] = clRmsXYZ[i] - clXYZ[i]*clXYZ[i]; + if(clRmsXYZ[i] > 0.0) { + clRmsXYZ[i] = TMath::Sqrt(clRmsXYZ[i]); + } else clRmsXYZ[i] = 0; + } else clRmsXYZ[i] = 0; + } + } else { + for(i=0; i<3; i++ ) { + clXYZ[i] = clRmsXYZ[i] = -1.; + } + } + // clRmsXYZ[i] ?? + if(phiSlope != 0.0 && logWeight > 0.0 && wtot) { + // Correction in phi direction (y - coords here); Aug 16; + // May be put to global level or seperate method + ycorr = clXYZ[1] * (1. + phiSlope); + //printf(" y %f : ycorr %f : slope %f \n", clXYZ[1], ycorr, phiSlope); + clXYZ[1] = ycorr; + } + fLocPos.SetX(clXYZ[0]); + fLocPos.SetY(clXYZ[1]); + fLocPos.SetZ(clXYZ[2]); + +// if (gDebug==2) +// printf("EvalLocalPosition: eta,phi,r = %f,%f,%f", fLocPos.X(), fLocPos.Y(), fLocPos.Z()) ; + fLocPosM = 0 ; // covariance matrix +} + +//_____________________________________________________________________________ +Bool_t AliEMCALRecPoint::EvalLocalPosition2(TClonesArray * digits, TArrayD &ed) +{ + // Evaluated local position of rec.point using digits + // and parametrisation of w0 and deff + //printf(" AliEMCALRecPoint::EvalLocalPosition2() \n"); + return AliEMCALRecPoint::EvalLocalPositionFromDigits(digits, ed, fLocPos); +} + +//_____________________________________________________________________________ +Bool_t AliEMCALRecPoint::EvalLocalPositionFromDigits(TClonesArray *digits, TArrayD &ed, TVector3 &locPos) +{ + // Used when digits should be recalibrated + static Double_t deff, w0, esum; + static Int_t iDigit; + // static AliEMCALDigit *digit; + + if(ed.GetSize() && (digits->GetEntries()!=ed.GetSize())) return kFALSE; + + // Calculate sum energy of digits + esum = 0.0; + for(iDigit=0; iDigitGetEntries(); iDigit++) { + digit = dynamic_cast(digits->At(iDigit)); + + geo->RelPosCellInSModule(digit->GetId(), idMax, deff, xyzi[0], xyzi[1], xyzi[2]); + + if(w0 > 0.0) w = TMath::Max( 0., w0 + TMath::Log(ed[iDigit] / esum)); + else w = ed[iDigit]; // just energy + + if(w>0.0) { + wtot += w ; + nstat++; + for(i=0; i<3; i++ ) { + clXYZ[i] += (w*xyzi[i]); + } + } + } + // cout << " wtot " << wtot << endl; + if (wtot > 0) { + for(i=0; i<3; i++ ) { + clXYZ[i] /= wtot; + } + locPos.SetX(clXYZ[0]); + locPos.SetY(clXYZ[1]); + locPos.SetZ(clXYZ[2]); + return kTRUE; + } else { + return kFALSE; + } + +} + +//_____________________________________________________________________________ +void AliEMCALRecPoint::GetDeffW0(const Double_t esum , Double_t &deff, Double_t &w0) +{ + // + // Aug 31, 2001 + // Applied for simulation data with threshold 3 adc + // Calculate efective distance (deff) and weigh parameter (w0) + // for coordinate calculation; 0.5 GeV < esum <100 GeV. + // Look to: http://rhic.physics.wayne.edu/~pavlinov/ALICE/SHISHKEBAB/RES/CALIB/GEOMCORR/deffandW0VaEgamma_2.gif + // + static Double_t e=0.0; + const Double_t kdp0=9.25147, kdp1=1.16700; // Hard coded now + const Double_t kwp0=4.83713, kwp1=-2.77970e-01, kwp2 = 4.41116; + + // No extrapolation here + e = esum<0.5?0.5:esum; + e = e>100.?100.:e; + + deff = kdp0 + kdp1*TMath::Log(e); + w0 = kwp0 / (1. + TMath::Exp(kwp1*(e+kwp2))); + //printf(" AliEMCALRecPoint::GetDeffW0 esum %5.2f : deff %5.2f : w0 %5.2f \n", esum, deff, w0); +} //______________________________________________________________________________ void AliEMCALRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits) { // This function calculates energy in the core, - // i.e. within a radius rad = 3cm around the center. Beyond this radius + // i.e. within a radius rad = fCoreEnergy around the center. Beyond this radius // in accordance with shower profile the energy deposition // should be less than 2% + // Unfinished - Nov 15,2006 + // Distance is calculate in (phi,eta) units AliEMCALDigit * digit ; - const Float_t kDeg2Rad = TMath::DegToRad() ; - // AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance(); Int_t iDigit; @@ -452,14 +795,16 @@ void AliEMCALRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits) EvalLocalPosition(logWeight, digits); } + Double_t phiPoint = fLocPos.Phi(), etaPoint = fLocPos.Eta(); + Double_t eta, phi, distance; for(iDigit=0; iDigit < fMulDigit; iDigit++) { digit = (AliEMCALDigit *) ( digits->At(fDigitsList[iDigit]) ) ; - Float_t etai = 0. ; - Float_t phii = 0. ; - geom->PosInAlice(digit->GetId(), etai, phii); - phii = phii * kDeg2Rad; + + eta = phi = 0.0; + fGeomPtr->EtaPhiFromIndex(digit->GetId(),eta, phi) ; + phi = phi * TMath::DegToRad(); - Float_t distance = TMath::Sqrt((etai-fLocPos.X())*(etai-fLocPos.X())+(phii-fLocPos.Y())*(phii-fLocPos.Y())) ; + distance = TMath::Sqrt((eta-etaPoint)*(eta-etaPoint)+(phi-phiPoint)*(phi-phiPoint)); if(distance < fCoreRadius) fCoreEnergy += fEnergyList[iDigit] ; } @@ -469,41 +814,39 @@ void AliEMCALRecPoint::EvalCoreEnergy(Float_t logWeight, TClonesArray * digits) void AliEMCALRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits) { // Calculates the axis of the shower ellipsoid in eta and phi + // in cell units - Double_t wtot = 0. ; + static TString gn(fGeomPtr->GetName()); + + Double_t wtot = 0.; Double_t x = 0.; Double_t z = 0.; Double_t dxx = 0.; Double_t dzz = 0.; Double_t dxz = 0.; - const Float_t kDeg2Rad = TMath::DegToRad(); - AliEMCALDigit * digit ; - - // AliEMCALGeometry * geom = AliEMCALGeometry::GetInstance(); - TString gn(geom->GetName()); + AliEMCALDigit * digit = 0; - Int_t iDigit; - - for(iDigit=0; iDigitAt(fDigitsList[iDigit]) ; - Float_t etai = 0. ; - Float_t phii = 0. ; - if(gn.Contains("SHISH")) { // have to be change - Feb 28, 2006 - //copied for shish-kebab geometry, ieta,iphi is cast as float as eta,phi conversion - // for this geometry does not exist - 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); - etai=(Float_t)ieta; - phii=(Float_t)iphi; - } else { - geom->EtaPhiFromIndex(digit->GetId(), etai, phii); - phii = phii * kDeg2Rad; - } + etai = phii = 0.; + // Nov 15,2006 - use cell numbers as coordinates + // Copied for shish-kebab geometry, ieta,iphi is cast as double as eta,phi + // We can use the eta,phi(or coordinates) of cell + nSupMod = nModule = nIphi = nIeta = iphi = ieta = 0; - Double_t w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; + fGeomPtr->GetCellIndex(digit->GetId(), nSupMod,nModule,nIphi,nIeta); + fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta); + etai=(Double_t)ieta; + phii=(Double_t)iphi; + + w = TMath::Max(0.,logWeight+TMath::Log(fEnergyList[iDigit]/fAmp ) ) ; + // fAmp summed amplitude of digits, i.e. energy of recpoint + // Gives smaller value of lambda than log weight + // w = fEnergyList[iDigit] / fAmp; // Nov 16, 2006 - try just energy dxx += w * etai * etai ; x += w * etai ; @@ -549,51 +892,54 @@ void AliEMCALRecPoint::EvalElipsAxis(Float_t logWeight,TClonesArray * digits) //______________________________________________________________________________ void AliEMCALRecPoint::EvalPrimaries(TClonesArray * digits) { - // Constructs the list of primary particles (tracks) which have contributed to this RecPoint + // Constructs the list of primary particles (tracks) which + // have contributed to this RecPoint and calculate deposited energy + // for each track AliEMCALDigit * digit ; - Int_t * tempo = new Int_t[fMaxTrack] ; + Int_t * primArray = new Int_t[fMaxTrack] ; + Float_t * dEPrimArray = new Float_t[fMaxTrack] ; Int_t index ; for ( index = 0 ; index < GetDigitsMultiplicity() ; index++ ) { // all digits digit = dynamic_cast(digits->At( fDigitsList[index] )) ; Int_t nprimaries = digit->GetNprimary() ; if ( nprimaries == 0 ) continue ; - Int_t * newprimaryarray = new Int_t[nprimaries] ; - Int_t ii ; - for ( ii = 0 ; ii < nprimaries ; ii++) - newprimaryarray[ii] = digit->GetPrimary(ii+1) ; - Int_t jndex ; for ( jndex = 0 ; jndex < nprimaries ; jndex++ ) { // all primaries in digit if ( fMulTrack > fMaxTrack ) { fMulTrack = fMaxTrack ; - Error("GetNprimaries", "increase fMaxTrack ") ; + Error("EvalPrimaries", "increase fMaxTrack ") ; break ; } - Int_t newprimary = newprimaryarray[jndex] ; + Int_t newPrimary = digit->GetPrimary(jndex+1); + Float_t dEPrimary = digit->GetDEPrimary(jndex+1); Int_t kndex ; Bool_t already = kFALSE ; for ( kndex = 0 ; kndex < fMulTrack ; kndex++ ) { //check if not already stored - if ( newprimary == tempo[kndex] ){ + if ( newPrimary == primArray[kndex] ){ already = kTRUE ; + dEPrimArray[kndex] += dEPrimary; break ; } } // end of check if ( !already && (fMulTrack < fMaxTrack)) { // store it - tempo[fMulTrack] = newprimary ; + primArray[fMulTrack] = newPrimary ; + dEPrimArray[fMulTrack] = dEPrimary ; fMulTrack++ ; } // store it } // all primaries in digit - delete [] newprimaryarray ; } // all digits - - fTracksList = new Int_t[fMulTrack] ; - for(index = 0; index < fMulTrack; index++) - fTracksList[index] = tempo[index] ; - - delete [] tempo ; + Int_t *sortIdx = new Int_t[fMulTrack]; + TMath::Sort(fMulTrack,dEPrimArray,sortIdx); + for(index = 0; index < fMulTrack; index++) { + fTracksList[index] = primArray[sortIdx[index]] ; + fDETracksList[index] = dEPrimArray[sortIdx[index]] ; + } + delete [] sortIdx; + delete [] primArray ; + delete [] dEPrimArray ; } @@ -603,60 +949,64 @@ void AliEMCALRecPoint::EvalParents(TClonesArray * digits) // Constructs the list of parent particles (tracks) which have contributed to this RecPoint AliEMCALDigit * digit ; - Int_t * tempo = new Int_t[fMaxParent] ; + Int_t * parentArray = new Int_t[fMaxTrack] ; + Float_t * dEParentArray = new Float_t[fMaxTrack] ; Int_t index ; for ( index = 0 ; index < GetDigitsMultiplicity() ; index++ ) { // all digits + if (fDigitsList[index] >= digits->GetEntries() || fDigitsList[index] < 0) + AliError(Form("Trying to get invalid digit %d (idx in WriteRecPoint %d)",fDigitsList[index],index)); digit = dynamic_cast(digits->At( fDigitsList[index] )) ; Int_t nparents = digit->GetNiparent() ; if ( nparents == 0 ) continue ; - Int_t * newparentarray = new Int_t[nparents] ; - Int_t ii ; - for ( ii = 0 ; ii < nparents ; ii++) - newparentarray[ii] = digit->GetIparent(ii+1) ; Int_t jndex ; for ( jndex = 0 ; jndex < nparents ; jndex++ ) { // all primaries in digit if ( fMulParent > fMaxParent ) { fMulTrack = - 1 ; - Error("GetNiparent", "increase fMaxParent") ; + Error("EvalParents", "increase fMaxParent") ; break ; } - Int_t newparent = newparentarray[jndex] ; + Int_t newParent = digit->GetIparent(jndex+1) ; + Float_t newdEParent = digit->GetDEParent(jndex+1) ; Int_t kndex ; Bool_t already = kFALSE ; for ( kndex = 0 ; kndex < fMulParent ; kndex++ ) { //check if not already stored - if ( newparent == tempo[kndex] ){ + if ( newParent == parentArray[kndex] ){ + dEParentArray[kndex] += newdEParent; already = kTRUE ; break ; } } // end of check - if ( !already && (fMulTrack < fMaxTrack)) { // store it - tempo[fMulParent] = newparent ; + if ( !already && (fMulParent < fMaxParent)) { // store it + parentArray[fMulParent] = newParent ; + dEParentArray[fMulParent] = newdEParent ; fMulParent++ ; } // store it } // all parents in digit - delete [] newparentarray ; } // all digits if (fMulParent>0) { - fParentsList = new Int_t[fMulParent] ; - for(index = 0; index < fMulParent; index++) - fParentsList[index] = tempo[index] ; + Int_t *sortIdx = new Int_t[fMulParent]; + TMath::Sort(fMulParent,dEParentArray,sortIdx); + for(index = 0; index < fMulParent; index++) { + fParentsList[index] = parentArray[sortIdx[index]] ; + fDEParentsList[index] = dEParentArray[sortIdx[index]] ; + } + delete [] sortIdx; } - delete [] tempo ; - + delete [] parentArray; + delete [] dEParentArray; } //____________________________________________________________________________ void AliEMCALRecPoint::GetLocalPosition(TVector3 & lpos) const { - // returns the position of the cluster in the local reference system of ALICE + // returns the position of the cluster in the local reference system + // of the sub-detector - lpos.SetX(fLocPos.X()) ; - lpos.SetY(fLocPos.Y()) ; - lpos.SetZ(fLocPos.Z()) ; + lpos = fLocPos; } //____________________________________________________________________________ @@ -664,9 +1014,56 @@ void AliEMCALRecPoint::GetGlobalPosition(TVector3 & gpos) const { // returns the position of the cluster in the global reference system of ALICE // These are now the Cartesian X, Y and Z - geom = AliEMCALGeometry::GetInstance(); // cout<<" geom "<GetGlobal(fLocPos, gpos, fSuperModuleNumber); + fGeomPtr->GetGlobal(fLocPos, gpos, fSuperModuleNumber); + +} + +//____________________________________________________________________________ +void AliEMCALRecPoint::GetGlobalPosition(TVector3 & gpos, TMatrixF & gmat) const +{ + // returns the position of the cluster in the global reference system of ALICE + // These are now the Cartesian X, Y and Z + // cout<<" geom "<GetGlobalEMCAL(this, gpos, gmat); + +} + +//_____________________________________________________________________________ +void AliEMCALRecPoint::EvalLocal2TrackingCSTransform() +{ + //Evaluates local to "tracking" c.s. transformation (B.P.). + //All evaluations should be completed before calling for this + //function. + //See ALICE PPR Chapter 5 p.18 for "tracking" c.s. definition, + //or just ask Jouri Belikov. :) + + SetVolumeId(AliGeomManager::LayerToVolUID(AliGeomManager::kEMCAL,GetSuperModuleNumber())); + + const TGeoHMatrix* tr2loc = GetTracking2LocalMatrix(); + if(!tr2loc) AliFatal(Form("No Tracking2LocalMatrix found.")); + + Double_t lxyz[3] = {fLocPos.X(),fLocPos.Y(),fLocPos.Z()}; + Double_t txyz[3] = {0,0,0}; + + tr2loc->MasterToLocal(lxyz,txyz); + SetX(txyz[0]); SetY(txyz[1]); SetZ(txyz[2]); + + if(AliLog::GetGlobalDebugLevel()>0) { + TVector3 gpos; TMatrixF gmat; + GetGlobalPosition(gpos,gmat); + Float_t gxyz[3]; + GetGlobalXYZ(gxyz); + AliInfo(Form("lCS-->(%.3f,%.3f,%.3f), tCS-->(%.3f,%.3f,%.3f), gCS-->(%.3f,%.3f,%.3f), gCScalc-\ +->(%.3f,%.3f,%.3f), supermodule %d", + fLocPos.X(),fLocPos.Y(),fLocPos.Z(), + GetX(),GetY(),GetZ(), + gpos.X(),gpos.Y(),gpos.Z(), + gxyz[0],gxyz[1],gxyz[2],GetSuperModuleNumber())); + } + } //____________________________________________________________________________ @@ -752,6 +1149,18 @@ Int_t AliEMCALRecPoint::GetNumberOfLocalMax(AliEMCALDigit ** maxAt, Float_t * } return iDigitN ; } + +//____________________________________________________________________________ +Int_t AliEMCALRecPoint::GetPrimaryIndex() const +{ + // Get the primary track index in TreeK which deposits the most energy + // in Digits which forms RecPoint. + + if (fMulTrack) + return fTracksList[0]; + return -12345; +} + //____________________________________________________________________________ void AliEMCALRecPoint::EvalTime(TClonesArray * digits){ // time is set to the time of the digit with the maximum energy @@ -790,6 +1199,26 @@ void AliEMCALRecPoint::Paint(Option_t *) gPad->PaintPolyMarker(1,&x,&y,"") ; } +//_____________________________________________________________________ +Double_t AliEMCALRecPoint::TmaxInCm(const Double_t e , const Int_t key) +{ + // e energy in GeV) + // key = 0(gamma, default) + // != 0(electron) + static Double_t ca = 4.82; // shower max parameter - first guess; ca=TMath::Log(1000./8.07) + static Double_t x0 = 1.23; // radiation lenght (cm) + static Double_t tmax = 0.; // position of electromagnetic shower max in cm + + tmax = 0.0; + if(e>0.1) { + tmax = TMath::Log(e) + ca; + if (key==0) tmax += 0.5; + else tmax -= 0.5; + tmax *= x0; // convert to cm + } + return tmax; +} + //______________________________________________________________________________ Float_t AliEMCALRecPoint::EtaToTheta(Float_t arg) const { @@ -805,10 +1234,10 @@ Float_t AliEMCALRecPoint::ThetaToEta(Float_t arg) const } //____________________________________________________________________________ -void AliEMCALRecPoint::Print(Option_t *) const +void AliEMCALRecPoint::Print(Option_t *opt) const { // Print the list of digits belonging to the cluster - return; + if(strlen(opt)==0) return; TString message ; message = "AliEMCALRecPoint:\n" ; message += " digits # = " ; @@ -835,11 +1264,22 @@ void AliEMCALRecPoint::Print(Option_t *) const printf("\n Local x %6.2f y %7.2f z %7.1f \n", fLocPos[0], fLocPos[1], fLocPos[2]); - message = " Multiplicity = %d" ; + message = " ClusterType = %d" ; + message += " Multiplicity = %d" ; message += " Cluster Energy = %f" ; message += " Core energy = %f" ; message += " Core radius = %f" ; message += " Number of primaries %d" ; message += " Stored at position %d" ; - Info("Print", message.Data(), fMulDigit, fAmp, fCoreEnergy, fCoreRadius, fMulTrack, GetIndexInList() ) ; + Info("Print", message.Data(), fClusterType, fMulDigit, fAmp, fCoreEnergy, fCoreRadius, fMulTrack, GetIndexInList() ) ; +} + +//___________________________________________________________ +Double_t AliEMCALRecPoint::GetPointEnergy() const +{ + //Returns energy .... + static double e; + e=0.0; + for(int ic=0; ic