Updating limits to cope both with p-A and A-p ZDC timing
[u/mrichter/AliRoot.git] / 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 }