X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=EMCAL%2FAliEMCALRecPoint.cxx;h=0b452642fc1f7f350c8ea51166eac5e299434988;hb=ea12701f2f35e3e60c95a672542f2de4ada91bc3;hp=96cc08165b69945efb70e6e4e3f33064554a140a;hpb=d8c2bd6979f96ecc6edf0c4ca3acce96ab4693e0;p=u%2Fmrichter%2FAliRoot.git diff --git a/EMCAL/AliEMCALRecPoint.cxx b/EMCAL/AliEMCALRecPoint.cxx index 96cc08165b6..0b452642fc1 100644 --- a/EMCAL/AliEMCALRecPoint.cxx +++ b/EMCAL/AliEMCALRecPoint.cxx @@ -59,7 +59,7 @@ AliEMCALRecPoint::AliEMCALRecPoint() fMaxDigit(100), fMulDigit(0), fMaxTrack(200), fMulTrack(0), fDigitsList(0), fTracksList(0), fClusterType(-1), fCoreEnergy(0), fDispersion(0), - fEnergyList(0), fTimeList(0), fAbsIdList(0), + fEnergyList(0), fAbsIdList(0), fTime(0.), fNExMax(0), fCoreRadius(10), //HG check this fDETracksList(0), fMulParent(0), fMaxParent(0), fParentsList(0), fDEParentsList(0), fSuperModuleNumber(0), @@ -81,7 +81,7 @@ AliEMCALRecPoint::AliEMCALRecPoint(const char *) 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]), + fEnergyList(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]), @@ -110,7 +110,7 @@ AliEMCALRecPoint::AliEMCALRecPoint(const AliEMCALRecPoint & rp) 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]), + fEnergyList(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]), @@ -124,8 +124,7 @@ AliEMCALRecPoint::AliEMCALRecPoint(const AliEMCALRecPoint & rp) for(Int_t i = 0; i < rp.fMulDigit; i++) { fEnergyList[i] = rp.fEnergyList[i]; - fTimeList[i] = rp.fTimeList[i]; - fAbsIdList[i] = rp.fAbsIdList[i]; + fAbsIdList[i] = rp.fAbsIdList[i]; } for(Int_t i = 0; i < rp.fMulTrack; i++) fDETracksList[i] = rp.fDETracksList[i]; @@ -142,8 +141,6 @@ AliEMCALRecPoint::~AliEMCALRecPoint() // dtor if ( fEnergyList ) delete[] fEnergyList ; - if ( fTimeList ) - delete[] fTimeList ; if ( fAbsIdList ) delete[] fAbsIdList ; if ( fDETracksList) @@ -164,37 +161,59 @@ AliEMCALRecPoint& AliEMCALRecPoint::operator= (const AliEMCALRecPoint &rp) if(&rp == this) return *this; - fGeomPtr = rp.fGeomPtr; - fAmp = rp.fAmp; + fGeomPtr = rp.fGeomPtr; + fAmp = rp.fAmp; fIndexInList = rp.fIndexInList; - fGlobPos = rp.fGlobPos; - fLocPos = rp.fLocPos; - fMaxDigit = rp.fMaxDigit; - fMulDigit = rp.fMulDigit; - fMaxTrack = rp.fMaxTrack; - fMulTrack = rp.fMaxTrack; + fGlobPos = rp.fGlobPos; + fLocPos = rp.fLocPos; + fMaxDigit = rp.fMaxDigit; + fMulDigit = rp.fMulDigit; + fMaxTrack = rp.fMaxTrack; + fMulTrack = rp.fMulTrack; + + if(fDigitsList) delete [] fDigitsList; + fDigitsList = new Int_t[rp.fMaxDigit]; + if(fTracksList) delete [] fTracksList; + fTracksList = new Int_t[rp.fMaxTrack]; for(Int_t i = 0; i= fMaxDigit ) { // increase the size of the lists fMaxDigit*=2 ; - Int_t * tempo = new Int_t[fMaxDigit]; + 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 * 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] ; + tempo [index] = fDigitsList[index] ; + tempoE [index] = fEnergyList[index] ; + tempoId[index] = fAbsIdList [index] ; } delete [] fDigitsList ; delete [] fEnergyList ; - delete [] fTimeList ; delete [] fAbsIdList ; fDigitsList = tempo; fEnergyList = tempoE; - fTimeList = tempoT; fAbsIdList = tempoId; } // if fDigitsList[fMulDigit] = digit.GetIndexInList() ; fEnergyList[fMulDigit] = energy ; - fTimeList[fMulDigit] = digit.GetTimeR() ; - fAbsIdList[fMulDigit] = digit.GetId(); + fAbsIdList [fMulDigit] = digit.GetId(); fMulDigit++ ; fAmp += energy ; if(shared) fSharedCluster = kTRUE; - - //GCB, May-2010, setting moved to EvalAll method, set the super module number for the largest energy digit position. - //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 @@ -288,15 +292,20 @@ Bool_t AliEMCALRecPoint::AreNeighbours(AliEMCALDigit * digit1, AliEMCALDigit * d // In case of a shared cluster, index of SM in C side, columns start at 48 and ends at 48*2-1 // C Side impair SM, nSupMod%2=1; A side pair SM nSupMod%2=0 if(fSharedCluster){ - if(nSupMod1%2) relid1[1]+=AliEMCALGeoParams::fgkEMCALCols; - else relid2[1]+=AliEMCALGeoParams::fgkEMCALCols; + //printf("Shared cluster in 2 SMs!\n"); + + // if(nSupMod1%2) relid1[1]+=AliEMCALGeoParams::fgkEMCALCols;//bad + // else relid2[1]+=AliEMCALGeoParams::fgkEMCALCols;//bad + if(nSupMod1%2) relid2[1]+=AliEMCALGeoParams::fgkEMCALCols; + else relid1[1]+=AliEMCALGeoParams::fgkEMCALCols; } rowdiff = TMath::Abs( relid1[0] - relid2[0] ) ; coldiff = TMath::Abs( relid1[1] - relid2[1] ) ; - if (( coldiff <= 1 ) && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0)) - areNeighbours = kTRUE ; + //if (( coldiff <= 1 ) && ( rowdiff <= 1 ) && (coldiff + rowdiff > 0)) + if ((coldiff + rowdiff == 1 )) + areNeighbours = kTRUE ; return areNeighbours; } @@ -330,28 +339,6 @@ Int_t AliEMCALRecPoint::Compare(const TObject * obj) const return rv ; } -// GCB, May-2010, Method not used, just comment it but remove? -//____________________________________________________________________________ -//Int_t AliEMCALRecPoint::DistancetoPrimitive(Int_t px, Int_t py) -//{ -// // Compute distance from point px,py to a AliEMCALRecPoint considered as a Tmarker -// // Compute the closest distance of approach from point px,py to this marker. -// // The distance is computed in pixels units. -// // HG Still need to update -> Not sure what this should achieve -// -// TVector3 pos(0.,0.,0.) ; -// GetLocalPosition(pos) ; -// Float_t x = pos.X() ; -// Float_t y = pos.Y() ; -// const Int_t kMaxDiff = 10; -// Int_t pxm = gPad->XtoAbsPixel(x); -// Int_t pym = gPad->YtoAbsPixel(y); -// Int_t dist = (px-pxm)*(px-pxm) + (py-pym)*(py-pym); -// -// if (dist > kMaxDiff) return 9999; -// return dist; -//} - //___________________________________________________________________________ void AliEMCALRecPoint::Draw(Option_t *option) { @@ -360,88 +347,8 @@ Int_t AliEMCALRecPoint::Compare(const TObject * obj) const AppendPad(option); } -// GCB, May-2010, Method not used, just comment it but remove? -//______________________________________________________________________________ -//void AliEMCALRecPoint::ExecuteEvent(Int_t /*event*/, Int_t, Int_t) -//{ -// // Execute action corresponding to one event -// // This member function is called when a AliEMCALRecPoint is clicked with the locator -// // -// // If Left button is clicked on AliEMCALRecPoint, the digits are switched on -// // and switched off when the mouse button is released. -// -// // static Int_t pxold, pyold; -// -// /* static TGraph * digitgraph = 0 ; -// static TPaveText* clustertext = 0 ; -// -// if (!gPad->IsEditable()) return; -// -// switch (event) { -// -// -// case kButton1Down:{ -// AliEMCALDigit * digit ; -// -// Int_t iDigit; -// Int_t relid[2] ; -// -// const Int_t kMulDigit=AliEMCALRecPoint::GetDigitsMultiplicity() ; -// Float_t * xi = new Float_t [kMulDigit] ; -// Float_t * zi = new Float_t [kMulDigit] ; -// -// for(iDigit = 0; iDigit < kMulDigit; iDigit++) { -// Fatal("AliEMCALRecPoint::ExecuteEvent", " -> Something wrong with the code"); -// digit = 0 ; //dynamic_cast((fDigitsList)[iDigit]); -// fGeomPtr->AbsToRelNumbering(digit->GetId(), relid) ; -// fGeomPtr->PosInAlice(relid, xi[iDigit], zi[iDigit]) ; -// } -// -// if (!digitgraph) { -// digitgraph = new TGraph(fMulDigit,xi,zi); -// digitgraph-> SetMarkerStyle(5) ; -// digitgraph-> SetMarkerSize(1.) ; -// digitgraph-> SetMarkerColor(1) ; -// digitgraph-> Draw("P") ; -// } -// if (!clustertext) { -// -// TVector3 pos(0.,0.,0.) ; -// GetLocalPosition(pos) ; -// clustertext = new TPaveText(pos.X()-10,pos.Z()+10,pos.X()+50,pos.Z()+35,"") ; -// Text_t line1[40] ; -// Text_t line2[40] ; -// sprintf(line1,"Energy=%1.2f GeV",GetEnergy()) ; -// sprintf(line2,"%d Digits",GetDigitsMultiplicity()) ; -// clustertext ->AddText(line1) ; -// clustertext ->AddText(line2) ; -// clustertext ->Draw(""); -// } -// gPad->Update() ; -// Print("") ; -// delete[] xi ; -// delete[] zi ; -// } -// -//break; -// -// case kButton1Up: -// if (digitgraph) { -// delete digitgraph ; -// digitgraph = 0 ; -// } -// if (clustertext) { -// delete clustertext ; -// clustertext = 0 ; -// } -// -// break; -// -// }*/ -//} - //____________________________________________________________________________ -void AliEMCALRecPoint::EvalAll(Float_t logWeight,TClonesArray * digits) +void AliEMCALRecPoint::EvalAll(Float_t logWeight,TClonesArray * digits, const Bool_t justClusters) { // Evaluates cluster parameters @@ -465,7 +372,10 @@ void AliEMCALRecPoint::EvalAll(Float_t logWeight,TClonesArray * digits) EvalParents(digits); //Called last because it sets the global position of the cluster? - EvalLocal2TrackingCSTransform(); + //Do not call it when recalculating clusters out of standard reconstruction + if(!justClusters){ + EvalLocal2TrackingCSTransform(); + } } @@ -544,196 +454,156 @@ void AliEMCALRecPoint::EvalDispersion(Float_t logWeight, TClonesArray * digits) //____________________________________________________________________________ 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())); + //For each EMC rec. point set the distance to the nearest bad channel. + //AliInfo(Form("%d bad channel(s) found.\n", caloped->GetDeadTowerCount())); //It is done in cell units and not in global or local position as before (Sept 2010) - - if(!caloped->GetDeadTowerCount()) return; - - //Get channels map of the supermodule where the cluster is. - TH2D* hMap = caloped->GetDeadMap(fSuperModuleNumber); - + + if(!caloped->GetDeadTowerCount()) return; + + //Get channels map of the supermodule where the cluster is. + TH2D* hMap = caloped->GetDeadMap(fSuperModuleNumber); + Int_t dRrow, dReta; - Float_t minDist = 10000.; - Float_t dist = 0.; + Float_t minDist = 10000.; + Float_t dist = 0.; Int_t nSupMod, nModule; Int_t nIphi, nIeta; Int_t iphi, ieta; fDigitIndMax = GetMaximalEnergyIndex(); fGeomPtr->GetCellIndex(fAbsIdList[fDigitIndMax], nSupMod,nModule,nIphi,nIeta); fGeomPtr->GetCellPhiEtaIndexInSModule(nSupMod,nModule,nIphi,nIeta, iphi,ieta); - - // 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; + 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); - + dRrow=TMath::Abs(irow-iphi); dReta=TMath::Abs(icol-ieta); dist=TMath::Sqrt(dRrow*dRrow+dReta*dReta); - if(dist < minDist) minDist = dist; + if(dist < minDist) minDist = dist; - // //Tower is bad, get the absId of the index. - // absId = fGeomPtr->GetAbsCellIdFromCellIndexes(fSuperModuleNumber, irow, icol); - // - // //Get the position of this tower. - // - // //Calculate the distance in local coordinates - // //fGeomPtr->RelPosCellInSModule(absId,cellpos); - // //Calculate distance between this tower and cluster, set if is smaller than previous. - // //dR = cellpos-fLocPos; - // - // //Calculate the distance in global coordinates - // fGeomPtr->GetGlobal(absId,cellpos); - // //Calculate distance between this tower and cluster, set if it is smaller than previous. - // dR = cellpos-fGlobPos; - // - // dist = dR.Mag(); - // if(dist < minDist) minDist = dist; - - } - } + } + } - //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module - if (fSharedCluster) { - TH2D* hMap2 = 0; - Int_t nSupMod2 = -1; - - //The only possible combinations are (0,1), (2,3) ... (10,11) - if(fSuperModuleNumber%2) nSupMod2 = fSuperModuleNumber-1; - else nSupMod2 = fSuperModuleNumber+1; - hMap2 = caloped->GetDeadMap(nSupMod2); + //In case the cluster is shared by 2 SuperModules, need to check the map of the second Super Module + if (fSharedCluster) { + TH2D* hMap2 = 0; + Int_t nSupMod2 = -1; - //Loop on tower status map of second super module - for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){ - for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){ - //Check if tower is bad. - if(hMap2->GetBinContent(icol,irow)==AliCaloCalibPedestal::kAlive) continue; - //printf("AliEMCALRecPoint::EvalDistanceToBadChannels() - Bad channel in SM %d, col %d, row %d\n",iSM,icol, irow); - + //The only possible combinations are (0,1), (2,3) ... (10,11) + if(fSuperModuleNumber%2) nSupMod2 = fSuperModuleNumber-1; + else nSupMod2 = fSuperModuleNumber+1; + hMap2 = caloped->GetDeadMap(nSupMod2); + + //Loop on tower status map of second super module + for(Int_t irow = 0; irow < AliEMCALGeoParams::fgkEMCALRows; irow++){ + for(Int_t icol = 0; icol < AliEMCALGeoParams::fgkEMCALCols; icol++){ + //Check if tower is bad. + if(hMap2->GetBinContent(icol,irow)==AliCaloCalibPedestal::kAlive) continue; + //printf("AliEMCALRecPoint::EvalDistanceToBadChannels() - Bad channel in SM %d, col %d, row %d\n",iSM,icol, irow); dRrow=TMath::Abs(irow-iphi); - + if(fSuperModuleNumber%2) { - dReta=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+ieta)); - } - else { - dReta=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-ieta); - } + dReta=TMath::Abs(icol-(AliEMCALGeoParams::fgkEMCALCols+ieta)); + } + else { + dReta=TMath::Abs(AliEMCALGeoParams::fgkEMCALCols+icol-ieta); + } - dist=TMath::Sqrt(dRrow*dRrow+dReta*dReta); + dist=TMath::Sqrt(dRrow*dRrow+dReta*dReta); if(dist < minDist) minDist = dist; - -// -// //Tower is bad, get the absId of the index. -// absId = fGeomPtr->GetAbsCellIdFromCellIndexes(nSupMod2, irow, icol); -// -// //Get the position of this tower. -// -// //Calculate the distance in global coordinates -// fGeomPtr->GetGlobal(absId,cellpos); -// //Calculate distance between this tower and cluster, set if it is smaller than previous. -// dR = cellpos-fGlobPos; -// -// dist = dR.Mag(); -// if(dist < minDist) minDist = dist; - } - } - - }// shared cluster in 2 SuperModules - - fDistToBadTower = minDist; - //printf("AliEMCALRecPoint::EvalDistanceToBadChannel() - Distance to Bad is %f cm, shared cluster? %d \n",fDistToBadTower,fSharedCluster); + + } + } + + }// shared cluster in 2 SuperModules + + fDistToBadTower = minDist; + //printf("AliEMCALRecPoint::EvalDistanceToBadChannel() - Distance to Bad is %f cm, shared cluster? %d \n",fDistToBadTower,fSharedCluster); } //____________________________________________________________________________ void AliEMCALRecPoint::EvalLocalPosition(Float_t logWeight, TClonesArray * digits) { - // Calculates the center of gravity in the local EMCAL-module coordinates - // Info("Print", " logWeight %f : cluster energy %f ", logWeight, fAmp); // for testing - - AliEMCALDigit * digit=0; - Int_t i=0, nstat=0; - - Double_t dist = TmaxInCm(Double_t(fAmp)); - //Int_t idMax = GetAbsIdMaxDigit();// idMax is not used at all in RelPosCellInSModule, why use it? - - Double_t clXYZ[3]={0.,0.,0.}, clRmsXYZ[3]={0.,0.,0.}, xyzi[3], wtot=0., w=0.; - - //printf(" dist : %f e : %f \n", dist, fAmp); - for(Int_t iDigit=0; iDigit(digits->At(fDigitsList[iDigit])) ; - + // Calculates the center of gravity in the local EMCAL-module coordinates + // Info("Print", " logWeight %f : cluster energy %f ", logWeight, fAmp); // for testing + + AliEMCALDigit * digit=0; + Int_t i=0, nstat=0; + + Double_t dist = TmaxInCm(Double_t(fAmp)); + //Int_t idMax = GetAbsIdMaxDigit();// idMax is not used at all in RelPosCellInSModule, why use it? + + Double_t clXYZ[3]={0.,0.,0.}, clRmsXYZ[3]={0.,0.,0.}, xyzi[3], wtot=0., w=0.; + + //printf(" dist : %f e : %f \n", dist, fAmp); + for(Int_t iDigit=0; iDigit(digits->At(fDigitsList[iDigit])) ; + if(!digit) { AliError("No Digit!!"); continue; } - //fGeomPtr->RelPosCellInSModule(digit->GetId(), idMax, dist, xyzi[0], xyzi[1], xyzi[2]); - fGeomPtr->RelPosCellInSModule(digit->GetId(), dist, xyzi[0], xyzi[1], xyzi[2]); - - //Temporal patch, due to mapping problem, need to swap "y" in one of the 2 SM, although no effect in position calculation. GCB 05/2010 - if(fSharedCluster && fSuperModuleNumber != fGeomPtr->GetSuperModuleNumber(digit->GetId())) xyzi[1]*=-1; - - //printf("EvalLocalPosition Cell: Id %i, SM %i : dist %f Local x,y,z %f %f %f \n", - // digit->GetId(), fGeomPtr->GetSuperModuleNumber(digit->GetId()), 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] ?? - -// // Cluster of one single digit, smear the position to avoid discrete position -// // smear x and z with +- 3 cm to uniform (avoid discrete effects). Tower size is approx 6 cm. -// // Rndm generates a number in ]0,1] -// if (fMulDigit==1) { -// clXYZ[0] += fGeomPtr->GetPhiTileSize()*(0.5 - gRandom->Rndm()); -// clXYZ[2] += fGeomPtr->GetEtaTileSize()*(0.5 - gRandom->Rndm()); -// } - - //Set position in local vector - fLocPos.SetX(clXYZ[0]); - fLocPos.SetY(clXYZ[1]); - fLocPos.SetZ(clXYZ[2]); - - if (gDebug==2) - printf("EvalLocalPosition Cluster: Local (x,y,z) = (%f,%f,%f) \n", fLocPos.X(), fLocPos.Y(), fLocPos.Z()) ; - + fGeomPtr->RelPosCellInSModule(digit->GetId(), dist, xyzi[0], xyzi[1], xyzi[2]); + + //Temporal patch, due to mapping problem, need to swap "y" in one of the 2 SM, although no effect in position calculation. GCB 05/2010 + if(fSharedCluster && fSuperModuleNumber != fGeomPtr->GetSuperModuleNumber(digit->GetId())) xyzi[1]*=-1; + + //printf("EvalLocalPosition Cell: Id %i, SM %i : dist %f Local x,y,z %f %f %f \n", + // digit->GetId(), fGeomPtr->GetSuperModuleNumber(digit->GetId()), 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.; + } + } + + // // Cluster of one single digit, smear the position to avoid discrete position + // // smear x and z with +- 3 cm to uniform (avoid discrete effects). Tower size is approx 6 cm. + // // Rndm generates a number in ]0,1] + // if (fMulDigit==1) { + // clXYZ[0] += fGeomPtr->GetPhiTileSize()*(0.5 - gRandom->Rndm()); + // clXYZ[2] += fGeomPtr->GetEtaTileSize()*(0.5 - gRandom->Rndm()); + // } + + //Set position in local vector + fLocPos.SetX(clXYZ[0]); + fLocPos.SetY(clXYZ[1]); + fLocPos.SetZ(clXYZ[2]); + + if (gDebug==2) + printf("EvalLocalPosition Cluster: Local (x,y,z) = (%f,%f,%f) \n", fLocPos.X(), fLocPos.Y(), fLocPos.Z()) ; + } @@ -761,7 +631,6 @@ void AliEMCALRecPoint::EvalGlobalPosition(Float_t logWeight, TClonesArray * digi } //Get the local coordinates of the cell - //fGeomPtr->RelPosCellInSModule(digit->GetId(), idMax, dist, lxyzi[0], lxyzi[1], lxyzi[2]); fGeomPtr->RelPosCellInSModule(digit->GetId(), dist, lxyzi[0], lxyzi[1], lxyzi[2]); //Now get the global coordinate @@ -800,7 +669,6 @@ void AliEMCALRecPoint::EvalGlobalPosition(Float_t logWeight, TClonesArray * digi clXYZ[i] = clRmsXYZ[i] = -1.; } } - // clRmsXYZ[i] ?? // // Cluster of one single digit, smear the position to avoid discrete position // // smear x and z with +- 3 cm to uniform (avoid discrete effects). Tower size is approx 6 cm. @@ -1110,7 +978,7 @@ void AliEMCALRecPoint::EvalPrimaries(TClonesArray * digits) // Constructs the list of primary particles (tracks) which // have contributed to this RecPoint and calculate deposited energy // for each track - + AliEMCALDigit * digit =0; Int_t * primArray = new Int_t[fMaxTrack] ; memset(primArray,-1,sizeof(Int_t)*fMaxTrack); @@ -1169,7 +1037,7 @@ void AliEMCALRecPoint::EvalPrimaries(TClonesArray * digits) void AliEMCALRecPoint::EvalParents(TClonesArray * digits) { // Constructs the list of parent particles (tracks) which have contributed to this RecPoint - + AliEMCALDigit * digit=0 ; Int_t * parentArray = new Int_t[fMaxTrack] ; memset(parentArray,-1,sizeof(Int_t)*fMaxTrack); @@ -1282,10 +1150,10 @@ void AliEMCALRecPoint::EvalLocal2TrackingCSTransform() SetX(txyz[0]); SetY(txyz[1]); SetZ(txyz[2]); if(AliLog::GetGlobalDebugLevel()>0) { - TVector3 gpos; //TMatrixF gmat; + TVector3 gpos; //TMatrixF gmat; //GetGlobalPosition(gpos,gmat); //Not doing anythin special, replace by next line. - fGeomPtr->GetGlobal(fLocPos, gpos, GetSuperModuleNumber()); - + fGeomPtr->GetGlobal(fLocPos, gpos, GetSuperModuleNumber()); + Float_t gxyz[3]; GetGlobalXYZ(gxyz); AliInfo(Form("lCS-->(%.3f,%.3f,%.3f), tCS-->(%.3f,%.3f,%.3f), gCS-->(%.3f,%.3f,%.3f), gCScalc-\ @@ -1371,6 +1239,7 @@ Int_t AliEMCALRecPoint::GetNumberOfLocalMax(AliEMCALDigit ** maxAt, Float_t * digit = maxAt[iDigit] ; for(iDigitN = 0; iDigitN < fMulDigit; iDigitN++) { + if(iDigitN == iDigit) continue;//the same digit digitN = (AliEMCALDigit *) digits->At(fDigitsList[iDigitN]) ; if ( AreNeighbours(digit, digitN) ) { @@ -1379,8 +1248,7 @@ Int_t AliEMCALRecPoint::GetNumberOfLocalMax(AliEMCALDigit ** maxAt, Float_t * // but may be digit too is not local max ? if(fEnergyList[iDigit] < fEnergyList[iDigitN] + locMaxCut) maxAt[iDigit] = 0 ; - } - else { + } else { maxAt[iDigit] = 0 ; // but may be digitN too is not local max ? if(fEnergyList[iDigit] > fEnergyList[iDigitN] - locMaxCut) @@ -1399,6 +1267,7 @@ Int_t AliEMCALRecPoint::GetNumberOfLocalMax(AliEMCALDigit ** maxAt, Float_t * iDigitN++ ; } } + return iDigitN ; }