PHOS/AliPHOSAodCluster.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 //  AliAODCaloCluster 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 "AliPHOSAodCluster.h"
  35 #include "AliPHOSCalibData.h"
  36 #include "AliAODCaloCells.h"
  37
  38 ClassImp(AliPHOSAodCluster)
  39
  40 //____________________________________________________________________________
  41 AliPHOSAodCluster::AliPHOSAodCluster() :
  42   AliAODCaloCluster(),fRecalibrated(0)
  43 {
  44   // ctor
  45 }
  46 //____________________________________________________________________________
  47 AliPHOSAodCluster::AliPHOSAodCluster(const AliAODCaloCluster & clu) :
  48   AliAODCaloCluster(clu),fRecalibrated(0)
  49 {
  50   // cpy ctor
  51 }
  52
  53 //____________________________________________________________________________
  54 AliPHOSAodCluster::~AliPHOSAodCluster()
  55 {
  56   // dtor
  57 }
  58 //____________________________________________________________________________
  59 void AliPHOSAodCluster::Recalibrate(AliPHOSCalibData * calibData,AliAODCaloCells *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   Double32_t * cellsAmpFraction = GetCellsAmplitudeFraction();
  74
  75   for(Int_t i=0; i < GetNCells(); i++){
  76     Int_t relId[4];
  77     phosgeom->AbsToRelNumbering(GetCellAbsId(i),relId) ;
  78     Int_t   module = relId[0];
  79     Int_t   column = relId[3];
  80     Int_t   row    = relId[2];
  81     Double_t energy = phsCells->GetCellAmplitude(GetCellAbsId(i)) ;
  82     cellsAmpFraction[i]*=energy*calibData->GetADCchannelEmc(module,column,row);
  83   }
  84
  85   SetCellsAmplitudeFraction(cellsAmpFraction);
  86   fRecalibrated=kTRUE;
  87 }
  88 //____________________________________________________________________________
  89 void  AliPHOSAodCluster::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 AliPHOSAodCluster::EvalEnergy(){
 100   if(!fRecalibrated) // no need to recalibrate
 101     return ;
 102
 103   Float_t energy=0. ;
 104   for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
 105     energy+=GetCellAmplitudeFraction(iDigit) ;
 106   }
 107
 108   SetE(energy);
 109   //Correct for nonlinearity later
 110
 111 }
 112 //____________________________________________________________________________
 113 void AliPHOSAodCluster::EnergyCorrection(){
 114   //apply nonlinearity correction same as in AliPHOSPIDv1.
 115   SetE(AliPHOSReconstructor::CorrectNonlinearity(E())) ;
 116 }
 117 //____________________________________________________________________________
 118 void AliPHOSAodCluster::EvalPID(AliPHOSPIDv1 * /*pid*/){
 119   //re-evaluate identification parameters
 120 //  pid->CalculatePID(E(),GetDispersion(),GetEmcCpvDistance(),GetTOF(),GetPID()) ;
 121 //  pid->CalculatePID(E(),GetDispersion(),GetM20(),GetM02(),GetEmcCpvDistance(),GetTOF(),GetPID()) ;
 122 }
 123 //____________________________________________________________________________
 124 void AliPHOSAodCluster::EvalCoord(Float_t logWeight, TVector3 &vtx)
 125 {
 126   // Calculates new center of gravity in the local PHOS-module coordinates
 127   // and tranfers into global ALICE coordinates
 128   // Calculates Dispersion and main axis
 129   if(!fRecalibrated) // no need to recalibrate
 130     return ;
 131
 132   Float_t wtot = 0. ;
 133   Int_t relid[4] ;
 134   Int_t phosMod=0 ;
 135   Float_t xMean = 0. ;
 136   Float_t zMean = 0. ;
 137
 138   AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance() ;
 139   if(!phosgeom)
 140     AliFatal("AliPHOSGeometry was not contructed\n") ;
 141
 142   for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
 143     Float_t xi ;
 144     Float_t zi ;
 145     phosgeom->AbsToRelNumbering(GetCellAbsId(iDigit), relid) ;
 146     phosgeom->RelPosInModule(relid, xi, zi);
 147     phosMod=relid[0] ;
 148     Double_t ei=GetCellAmplitudeFraction(iDigit) ;
 149     if (E() > 0 && ei > 0) {
 150       Float_t w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ;
 151       xMean+= xi * w ;
 152       zMean+= zi * w ;
 153       wtot += w ;
 154     }
 155     else
 156       AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E()));
 157   }
 158   if (wtot>0) {
 159     xMean /= wtot ;
 160     zMean /= wtot ;
 161   }
 162   else
 163     AliError(Form("Wrong weight %f\n", wtot));
 164
 165
 166 // Calculates the dispersion and second momenta
 167   Double_t d=0. ;
 168   Double_t dxx  = 0.;
 169   Double_t dzz  = 0.;
 170   Double_t dxz  = 0.;
 171   for(Int_t iDigit=0; iDigit < GetNCells(); iDigit++) {
 172     Float_t xi ;
 173     Float_t zi ;
 174     phosgeom->AbsToRelNumbering(GetCellAbsId(iDigit), relid) ;
 175     phosgeom->RelPosInModule(relid, xi, zi);
 176     Double_t ei=GetCellAmplitudeFraction(iDigit) ;
 177     if (E() > 0 && ei > 0) {
 178       Float_t w = TMath::Max( 0., logWeight + TMath::Log(ei/E()) ) ;
 179       d += w*((xi-xMean)*(xi-xMean) + (zi-zMean)*(zi-zMean) ) ;
 180       dxx  += w * xi * xi ;
 181       dzz  += w * zi * zi ;
 182       dxz  += w * xi * zi ;
 183    }
 184     else
 185       AliError(Form("Wrong energy %f and/or amplitude %f\n", ei, E()));
 186   }
 187
 188   if (wtot>0) {
 189     d /= wtot ;
 190     dxx /= wtot ;
 191     dzz /= wtot ;
 192     dxz /= wtot ;
 193     dxx -= xMean * xMean ;
 194     dzz -= zMean * zMean ;
 195     dxz -= xMean * zMean ;
 196     SetM02(0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )) ;
 197     SetM20(0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz ));
 198   }
 199   else{
 200     AliError(Form("Wrong weight %f\n", wtot));
 201     d=0. ;
 202     SetM20(0.) ;
 203     SetM02(0.) ;
 204   }
 205
 206   if (d>=0)
 207     SetDispersion(TMath::Sqrt(d)) ;
 208   else
 209     SetDispersion(0) ;
 210
 211
 212   // Correction for the depth of the shower starting point (TDR p 127)
 213   Float_t para = 0.925 ;
 214   Float_t parb = 6.52 ;
 215
 216   TVector3 vInc ;
 217   phosgeom->GetIncidentVector(vtx,phosMod,xMean,zMean,vInc) ;
 218
 219   Float_t depthx = 0.;
 220   Float_t depthz = 0.;
 221   if (E()>0&&vInc.Y()!=0.) {
 222     depthx = ( para * TMath::Log(E()) + parb ) * vInc.X()/TMath::Abs(vInc.Y()) ;
 223     depthz = ( para * TMath::Log(E()) + parb ) * vInc.Z()/TMath::Abs(vInc.Y()) ;
 224   }
 225   else
 226     AliError(Form("Wrong amplitude %f\n", E()));
 227
 228   xMean-= depthx  ;
 229   zMean-= depthz  ;
 230
 231   //Go to the global system
 232   TVector3 gps ;
 233   phosgeom->Local2Global(phosMod, xMean, zMean, gps) ;
 234   SetPositionAt(gps[0],0) ;
 235   SetPositionAt(gps[1],1) ;
 236   SetPositionAt(gps[2],2) ;
 237 }