+//____________________________________________________________________________
+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; iDigit<digits->GetEntries(); iDigit++) {
+ digit = dynamic_cast<AliEMCALDigit *>(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(" <I> 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; iDigit<ed.GetSize(); iDigit++) esum += ed[iDigit];
+
+ GetDeffW0(esum, deff, w0);
+
+ return EvalLocalPositionFromDigits(esum, deff, w0, digits, ed, locPos);
+}
+
+Bool_t AliEMCALRecPoint::EvalLocalPositionFromDigits(const Double_t esum, const Double_t deff, const Double_t w0,
+ TClonesArray *digits, TArrayD &ed, TVector3 &locPos)
+{
+ static AliEMCALDigit *digit;
+
+ Int_t i=0, nstat=0, idMax=-1;
+ Double_t clXYZ[3]={0.,0.,0.}, xyzi[3], wtot=0., w=0.;
+
+ AliEMCALGeometry* geo = AliEMCALGeometry::GetInstance(); // Get pointer to EMCAL geometry
+
+ for(Int_t iDigit=0; iDigit<digits->GetEntries(); iDigit++) {
+ digit = dynamic_cast<AliEMCALDigit *>(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 dp0=9.25147, dp1=1.16700; // Hard coded now
+ const Double_t wp0=4.83713, wp1=-2.77970e-01, wp2 = 4.41116;
+
+ // No extrapolation here
+ e = esum<0.5?0.5:esum;
+ e = e>100.?100.:e;
+
+ deff = dp0 + dp1*TMath::Log(e);
+ w0 = wp0 / (1. + TMath::Exp(wp1*(e+wp2)));
+ //printf("<I> AliEMCALRecPoint::GetDeffW0 esum %5.2f : deff %5.2f : w0 %5.2f \n", esum, deff, w0);
+}