PHOS/AliPHOSEsdCluster.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //_________________________________________________________________________
  17 //  AliESDCaloCluster extension for PHOS to recalculate cluster
  18 //  parameters in case of recalibration.
  19 //*--
  20 //*-- Author: Dmitri Peressounko (RRC KI)
  21
  22
  23 // --- ROOT system ---
  24 #include "TVector3.h"
  25 #include "TMath.h"
  26
  27 // --- Standard library ---
  28
  29 // --- AliRoot header files ---
  30 #include "AliLog.h"
  31 #include "AliPHOSGeometry.h"
  32 #include "AliPHOSPIDv1.h"
  33 #include "AliPHOSReconstructor.h"
  34 #include "AliPHOSEsdCluster.h"
  35 #include "AliPHOSCalibData.h"
  36 #include "AliESDCaloCells.h"
  37
  38 ClassImp(AliPHOSEsdCluster)
  39
  40 //____________________________________________________________________________
  41 AliPHOSEsdCluster::AliPHOSEsdCluster() :
  42   AliESDCaloCluster(),fRecalibrated(0)
  43 {
  44   // ctor
  45 }
  46 //____________________________________________________________________________
  47 AliPHOSEsdCluster::AliPHOSEsdCluster(const AliESDCaloCluster & clu) :
  48   AliESDCaloCluster(clu),fRecalibrated(0)
  49 {
  50   // cpy ctor
  51 }
  52
  53 //____________________________________________________________________________
  54 AliPHOSEsdCluster::~AliPHOSEsdCluster()
  55 {
  56   // dtor
  57 }
  58 //____________________________________________________________________________
  59 void AliPHOSEsdCluster::Recalibrate(AliPHOSCalibData * calibData,AliESDCaloCells *phsCells){
  60   //If not done yet, apply recalibration coefficients to energies list
  61   //NOTE that after recalibration fCellsAmpFraction contains not FRACTION but FULL energy
  62
  63   if(fRecalibrated)
  64    return ;
  65
  66   if(!calibData)
  67     return ;
  68
  69   AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance() ;
  70   if(!phosgeom)
  71     AliFatal("AliPHOSGeometry was not contructed\n") ;
  72
  73   for(Int_t i=0; i<fNCells; i++){
  74     Int_t relId[4];
  75     phosgeom->AbsToRelNumbering(fCellsAbsId[i],relId) ;
  76     Int_t   module = relId[0];
  77     Int_t   column = relId[3];
  78     Int_t   row    = relId[2];
  79     Double_t energy = phsCells->GetCellAmplitude(fCellsAbsId[i]) ;
  80 Double_t fr = fCellsAmpFraction[i] ;
  81     fCellsAmpFraction[i]*=energy*calibData->GetADCchannelEmc(module,column,row);
  82     printf("      mod=%d, col=%d, row=%d \n",module,column,row);
  83     printf("  cell=%d, Amp=%f, E(i)=%f, Fr=%f, ci=%f \n",i,energy,fCellsAmpFraction[i],fr,calibData->GetADCchannelEmc(module,column,row)) ;
  84   }
  85 printf("----\n") ;
  86   fRecalibrated=kTRUE;
  87 }
  88 //____________________________________________________________________________
  89 void  AliPHOSEsdCluster::EvalAll(Float_t logWeight, TVector3 &vtx){
  90     //If recalibrated - recalculate all cluster parameters
  91   if(!fRecalibrated)
  92     return ;
  93
  94   EvalEnergy() ; //Energy should be evaluated first
  95   EvalCoord(logWeight, vtx) ;
  96
  97 }
  98 //____________________________________________________________________________
  99 void AliPHOSEsdCluster::EvalEnergy(){
 100   if(!fRecalibrated) // no need to recalibrate
 101     return ;
 102
 103   fEnergy=0. ;
 104   for(Int_t iDigit=0; iDigit<fNCells; iDigit++) {
 105     fEnergy+=fCellsAmpFraction[iDigit] ;
 106   }
 107   //Correct for nonlinearity later
 108   if(fEnergy==0.){
 109     printf("fEnergy=0 \n") ;
 110     printf("n Digits = %d \n",fNCells) ;
 111     for(Int_t iDigit=0; iDigit<fNCells; iDigit++) {
 112       printf("E(%d)=%f \n",iDigit,fCellsAmpFraction[iDigit]);
 113     }
 114   }
 115
 116 }
 117 //____________________________________________________________________________
 118 void AliPHOSEsdCluster::EnergyCorrection(){
 119   //apply nonlinearity correction
 120   fEnergy = AliPHOSReconstructor::CorrectNonlinearity(fEnergy) ;
 121 }
 122 //____________________________________________________________________________
 123 void AliPHOSEsdCluster::EvalPID(AliPHOSPIDv1 * /*pid*/){
 124   //re-evaluate identification parameters
 125 //  pid->CalculatePID(fEnergy,fDispersion,fEmcCpvDistance,tof,fPID) ;
 126 //  pid->CalculatePID(fEnergy,fDispersion,fM20,fM02,fEmcCpvDistance,tof,fPID) ;
 127 }
 128 //____________________________________________________________________________
 129 void AliPHOSEsdCluster::EvalCoord(Float_t logWeight, TVector3 &vtx)
 130 {
 131   // Calculates new center of gravity in the local PHOS-module coordinates
 132   // and tranfers into global ALICE coordinates
 133   // Calculates Dispersion and main axis
 134   if(!fRecalibrated) // no need to recalibrate
 135     return ;
 136
 137   Float_t wtot = 0. ;
 138   Int_t relid[4] ;
 139   Int_t phosMod=0 ;
 140   Float_t xMean = 0. ;
 141   Float_t zMean = 0. ;
 142
 143   AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance() ;
 144   if(!phosgeom)
 145     AliFatal("AliPHOSGeometry was not contructed\n") ;
 146
 147   for(Int_t iDigit=0; iDigit<fNCells; iDigit++) {
 148     Float_t xi ;
 149     Float_t zi ;
 150     phosgeom->AbsToRelNumbering(fCellsAbsId[iDigit], relid) ;
 151     phosgeom->RelPosInModule(relid, xi, zi);
 152     phosMod=relid[0] ;
 153     Double_t ei=fCellsAmpFraction[iDigit] ;
 154     if (fEnergy>0 && ei>0) {
 155       Float_t w = TMath::Max( 0., logWeight + TMath::Log(ei/fEnergy) ) ;
 156       xMean+= xi * w ;
 157       zMean+= zi * w ;
 158       wtot += w ;
 159     }
 160     else
 161       AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, fEnergy));
 162   }
 163   if (wtot>0) {
 164     xMean /= wtot ;
 165     zMean /= wtot ;
 166   }
 167   else
 168     AliError(Form("Wrong weight %f\n", wtot));
 169
 170
 171 // Calculates the dispersion and second momenta
 172   Double_t d=0. ;
 173   Double_t dxx  = 0.;
 174   Double_t dzz  = 0.;
 175   Double_t dxz  = 0.;
 176   for(Int_t iDigit=0; iDigit < fNCells; iDigit++) {
 177     Float_t xi ;
 178     Float_t zi ;
 179     phosgeom->AbsToRelNumbering(fCellsAbsId[iDigit], relid) ;
 180     phosgeom->RelPosInModule(relid, xi, zi);
 181     Double_t ei=fCellsAmpFraction[iDigit] ;
 182     if (fEnergy>0 && ei>0) {
 183       Float_t w = TMath::Max( 0., logWeight + TMath::Log(ei/fEnergy) ) ;
 184       d += w*((xi-xMean)*(xi-xMean) + (zi-zMean)*(zi-zMean) ) ;
 185       dxx  += w * xi * xi ;
 186       dzz  += w * zi * zi ;
 187       dxz  += w * xi * zi ;
 188    }
 189     else
 190       AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, fEnergy));
 191   }
 192
 193   if (wtot>0) {
 194     d /= wtot ;
 195     dxx /= wtot ;
 196     dzz /= wtot ;
 197     dxz /= wtot ;
 198     dxx -= xMean * xMean ;
 199     dzz -= zMean * zMean ;
 200     dxz -= xMean * zMean ;
 201     fM02 =  0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
 202     fM20 =  0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
 203   }
 204   else{
 205     AliError(Form("Wrong weight %f\n", wtot));
 206     d=0. ;
 207     fM20=0. ;
 208     fM02=0. ;
 209   }
 210
 211   if (d>=0)
 212     fDispersion = TMath::Sqrt(d) ;
 213   else
 214     fDispersion = 0 ;
 215
 216
 217   // Correction for the depth of the shower starting point (TDR p 127)
 218   Float_t para = 0.925 ;
 219   Float_t parb = 6.52 ;
 220
 221   TVector3 vInc ;
 222   phosgeom->GetIncidentVector(vtx,phosMod,xMean,zMean,vInc) ;
 223
 224   Float_t depthx = 0.;
 225   Float_t depthz = 0.;
 226   if (fEnergy>0&&vInc.Y()!=0.) {
 227     depthx = ( para * TMath::Log(fEnergy) + parb ) * vInc.X()/TMath::Abs(vInc.Y()) ;
 228     depthz = ( para * TMath::Log(fEnergy) + parb ) * vInc.Z()/TMath::Abs(vInc.Y()) ;
 229   }
 230   else
 231     AliError(Form("Wrong amplitude %f\n", fEnergy));
 232
 233   xMean-= depthx  ;
 234   zMean-= depthz  ;
 235
 236   //Go to the global system
 237   TVector3 gps ;
 238   phosgeom->Local2Global(phosMod, xMean, zMean, gps) ;
 239   fGlobalPos[0]=gps[0] ;
 240   fGlobalPos[1]=gps[1] ;
 241   fGlobalPos[2]=gps[2] ;
 242 }