Protections for coverity: DIVIDE_BY_ZERO
[u/mrichter/AliRoot.git] / JETAN / AliJetHadronCorrectionv1.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2002, 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 // To be modified for hadron correction using particle ID and track merging
18 // Author : magali.estienne@subatech.in2p3.fr
19 //===============================================================
20 // Author : Mark Horner (LBL/UCT)
21
22 // --- Standard library ---
23 #include "Riostream.h"
24 #include "TMath.h"
25
26 // --- AliRoot header files ---
27 #include "AliAODTrack.h"
28 #include "AliJetDummyGeo.h"
29 #include "AliJetHadronCorrectionv1.h"
30
31 static Double_t etaGrid[HCPARAMETERS]={ 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.67};
32
33 ClassImp(AliJetHadronCorrectionv1)
34
35 Double_t AliJetHadronCorrectionv1::fgParLookup[HCPARAMETERS][HCPARAMETERSETS] = 
36 {  
37
38     {7.52624e-2  , 9.80119e-2   , 1.29086e-2},
39     {-2.07449e-2 , 2.22951e-2   , 8.31971e-2},
40     {-3.64861e-4 , -3.65680e-03 , -3.07148e-02},
41     {5.35252e-3  , 5.12234e-03  , 2.25559e-02},
42     {1.27106e-1  , 1.25248e-01  , 1.06352e-01},
43     {-6.72909e-2 , 9.78677e-01  , -1.35909e-01},
44     {-3.80660e-5 , -2.77522e-05 , -2.05218e-05},
45     {1.54855e-7  , 1.00604e-07  , 1.02481e-07}
46 };
47
48 AliJetHadronCorrectionv1* AliJetHadronCorrectionv1::fgHadrCorr = 0;
49
50 AliJetHadronCorrectionv1::AliJetHadronCorrectionv1(const char *name,const char *title) 
51     :AliJetHadronCorrection(name, title), 
52      fSamplingFraction(0)
53 {
54   fgHadrCorr = this;
55   for (Int_t i = 0; i < 8; i++) fPar[i] = 0.;
56 }
57
58 AliJetHadronCorrectionv1*
59 AliJetHadronCorrectionv1::Instance()
60 {
61   // return pointer to global instance. Instantiate if needed
62   if (! fgHadrCorr) fgHadrCorr = new AliJetHadronCorrectionv1();
63   return fgHadrCorr;
64 }
65
66 void AliJetHadronCorrectionv1::SetGeometry2(const AliJetDummyGeo *geometry)
67 {
68   // Initialise EMCAL geometry
69     if (!geometry)
70     {
71         SetParameters();
72 //      commented out since geometry is 0x0 after null check
73 //      fSamplingFraction = geometry->GetSampling();
74     }else
75     {
76         SetParameters(geometry->GetName());
77         fSamplingFraction = geometry->GetSampling();
78     }
79     return;     
80 }       
81
82 /*
83 ###########################################################################
84 ###########################################################################
85
86 This will have to be modified with the inclusion of the new EMCAL geometry
87 That means I guess a new study of the hadron correction...
88
89 Geometry : SHISH...
90
91 ###########################################################################
92 ###########################################################################
93 */
94
95 void AliJetHadronCorrectionv1::SetGeometry(TString /*name*/,Double_t fs)
96 {
97   // Initialise EMCAL geometry
98   fSamplingFraction = fs;       
99
100   /*
101   if ( name == ""              ||
102        name == "EMCAL_5655_21" ||
103        name == "EMCALArch1a"   ||
104        name == "EMCALArch1aN"  ||
105        name == "G56_2_55_19"   ||
106        name == "G56_2_55_19_104_14" )
107     { // set parameters to this hadron correction
108       cout<<"HC parameters!"<<endl;
109       for (Int_t i=0;i<6;i++)
110         {
111           fPar[i] = fgParLookup[i][0];  
112           cout <<fPar[i]<<endl;
113         }
114     }else if( name == "EMCAL_6564_21" ||  
115               name == "G65_2_64_19" )
116     {     
117       cout<<"HC parameters!"<<endl;
118       for (Int_t i=0;i<6;i++)
119         {
120           fPar[i] = fgParLookup[i][1];  
121           cout <<fPar[i]<<endl;
122         }
123     }else
124     {
125       printf("Geometry not defined in hadron correction\n"); 
126     }     
127   */
128
129 }       
130
131 // Double_t AliJetHadronCorrectionv1::GetEnergy(Double_t pmom, Double_t eta, Int_t /*gid*/)
132 Double_t AliJetHadronCorrectionv1::GetEnergy(Double_t pmom, Double_t eta, Int_t /*gid*/)
133 {
134   // Return parametrised energy response
135   Double_t etai = TMath::Abs(eta); 
136
137   if(etai < etaGrid[1]) {
138     for (Int_t i=0;i<8;i++)
139       {
140         fPar[i] = fgParLookup[i][1];  
141         //      cout << "fPar[" << i << "]: " << fPar[i] << endl;
142       }
143   } else if(etai >= etaGrid[1] && etai <= etaGrid[HCPARAMETERS-2]) {
144     for (Int_t i=0;i<8;i++)
145       {
146         fPar[i] = fgParLookup[i][0];  
147         //      cout << "fPar[" << i << "]: " << fPar[i] << endl;
148       }
149   } else {
150     for (Int_t i=0;i<8;i++)
151       {
152         fPar[i] = fgParLookup[i][2];  
153         //      cout << "fPar[" << i << "]: " << fPar[i] << endl;
154       }
155
156   }
157
158   Double_t value = fPar[5]*pow(etai,3) + 
159                    fPar[0]*pow(etai,2) + 
160                    fPar[1]*etai        +
161                    fPar[2]*etai*pmom   +
162                    fPar[3]*pmom        + 
163                    fPar[4]             +
164                    fPar[6]*pow(pmom,2) +
165                    fPar[7]*pow(pmom,3);
166
167   return fSamplingFraction*value;
168    
169 }
170
171 void AliJetHadronCorrectionv1::TrackPositionEMCal(const AliAODTrack* track,Double_t &eta, Double_t &phi)
172 {
173 // Return track position on EMCal
174   AliAODPid*    pid   = (AliAODPid*) track->GetDetPid();
175     
176   if(pid) {
177     Double_t emcpos[3];
178     pid->GetEMCALPosition(emcpos);      
179     TVector3 tpos(emcpos[0],emcpos[1],emcpos[2]);
180     
181     eta = tpos.Eta();
182     phi = tpos.Phi();
183
184   }
185
186 }