Coverity fix
[u/mrichter/AliRoot.git] / TRD / TRDbase / AliTRDdEdxReconUtils.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 // TRD dEdx recon utils
18 // xx
19 // xx
20 // xx
21 // xx
22 //
23 //  Xianguo Lu 
24 //  lu@physi.uni-heidelberg.de
25 //  Xianguo.Lu@cern.ch
26 //  
27 //
28
29 #include "TF1.h"
30 #include "TFile.h"
31 #include "TH1D.h"
32 #include "TH2D.h"
33 #include "THnSparse.h"
34 #include "TMath.h"
35 #include "TMatrixD.h"
36 #include "TMinuit.h"
37 #include "TObjArray.h"
38 #include "TRandom3.h"
39 #include "TStopwatch.h"
40 #include "TVectorD.h"
41
42 #include "TTreeStream.h"
43
44 #include "AliCDBId.h"
45 #include "AliCDBMetaData.h"
46 #include "AliCDBStorage.h"
47 #include "AliESDEvent.h"
48 #include "AliESDfriendTrack.h"
49 #include "AliESDtrack.h"
50 #include "AliTRDcalibDB.h"
51 #include "AliTRDCalROC.h"
52 #include "AliTRDtrackV1.h"
53
54 #include "AliTRDdEdxBaseUtils.h"
55 #include "AliTRDdEdxReconUtils.h"
56
57 #define EPSILON 1e-8
58
59 Int_t AliTRDdEdxReconUtils::ApplyCalib(const Int_t nc0, TVectorD *arrayQ, TVectorD *arrayX, const TObjArray *cobj)
60 {
61   //
62   //apply calibration on arrayQ
63   //
64   if(!cobj){ printf("AliTRDdEdxReconUtils::ApplyCalib error gain array null!!\n"); exit(1);}
65
66   TVectorD tmpq(arrayQ->GetNrows());
67   TVectorD tmpx(arrayX->GetNrows());
68   Int_t ncls = 0;
69
70   const TVectorD * gain = (TVectorD*) cobj->At(0); 
71   for(Int_t ii=0; ii<nc0; ii++){
72     const Double_t dq = (*arrayQ)[ii];
73     const Int_t xx = (Int_t)(*arrayX)[ii];
74     const Double_t gg = (*gain)[xx];
75
76     if(gg<EPSILON){
77       continue;
78     }
79
80     tmpq[ncls] = dq*gg;
81     tmpx[ncls] = xx;
82     ncls++;
83   }
84
85   (*arrayQ)=tmpq;
86   (*arrayX)=tmpx;
87
88   return ncls;
89 }
90
91 Double_t AliTRDdEdxReconUtils::ToyCook(const Bool_t kinvq, Int_t &ncluster, TVectorD *arrayQ, TVectorD *arrayX, const TObjArray *cobj)
92 {
93   //
94   //template for cookdedx
95   //
96   if(cobj){
97     if(arrayQ && arrayX){
98       ncluster = ApplyCalib(ncluster, arrayQ, arrayX, cobj);
99     }
100     else{
101       printf("AliTRDdEdxReconUtils::ToyCook arrayQ arrayX null, applycalib can not be applied!\n"); exit(1);
102     }
103   }
104
105   Double_t lowFrac =-999, highFrac = -999;
106   if(kinvq){
107     lowFrac = AliTRDdEdxBaseUtils::Q1Frac(); highFrac = 0.99;
108   }
109   else{
110     lowFrac = 0.01; highFrac = AliTRDdEdxBaseUtils::Q0Frac();
111   }
112
113   Double_t meanQ = AliTRDdEdxBaseUtils::TruncatedMean(ncluster, arrayQ->GetMatrixArray(), lowFrac, highFrac);
114   if(kinvq){
115     if(meanQ>EPSILON){
116       meanQ = 1/meanQ;
117     }
118   }
119
120   return meanQ;
121 }
122
123 Double_t AliTRDdEdxReconUtils::CombineddEdx(const Bool_t kinvq, Int_t &concls, TVectorD *coarrayQ, TVectorD *coarrayX, const Int_t tpcncls, const TVectorD *tpcarrayQ, const TVectorD *tpcarrayX, const Int_t trdncls, const TVectorD *trdarrayQ, const TVectorD *trdarrayX)
124 {
125   //
126   //combine tpc and trd dedx
127   //
128
129   for(Int_t iq=0; iq<tpcncls; iq++){
130     (*coarrayQ)[iq]=(*tpcarrayQ)[iq];
131     if(tpcarrayX && trdarrayX && coarrayX){
132       (*coarrayX)[iq]=(*tpcarrayX)[iq];
133     }
134   }
135   for(Int_t iq=0; iq<trdncls; iq++){
136     (*coarrayQ)[tpcncls+iq]=(*trdarrayQ)[iq];
137     if(tpcarrayX && trdarrayX && coarrayX){
138       (*coarrayX)[tpcncls+iq]=159+(*trdarrayX)[iq];
139     }
140   }
141
142   concls=trdncls+tpcncls;
143
144   const Double_t coQ = ToyCook(kinvq, concls, coarrayQ, coarrayX);
145
146   return coQ;
147 }
148
149 Double_t AliTRDdEdxReconUtils::GetPadGain(const Int_t det, const Int_t icol, const Int_t irow)
150 {
151   //
152   //get pad calibration
153   //
154   AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
155   if (!calibration) {
156     printf("AliTRDdEdxReconUtils::GetPadCalib No AliTRDcalibDB instance available\n"); exit(1);
157   }
158   AliTRDCalROC * calGainFactorROC = calibration->GetGainFactorROC(det);
159   if(!calGainFactorROC){
160     printf("AliTRDdEdxReconUtils::GetPadCalib no calGainFactorROC!\n"); exit(1);
161   }
162
163   Double_t padgain = -999;
164   if( icol >= 0 && 
165       icol < calGainFactorROC->GetNcols() && 
166       irow >=0 && 
167       irow < calGainFactorROC->GetNrows()){
168     padgain = calGainFactorROC->GetValue(icol, irow);
169     if(padgain<EPSILON){
170       printf("AliTRDdEdxReconUtils::GetPadGain padgain 0! %f %f -- %d %d %d -- %d %d\n", padgain, EPSILON, det, icol, irow, calGainFactorROC->GetNcols(), calGainFactorROC->GetNrows()); exit(1);
171     }
172   }
173   else{
174     //printf("\nAliTRDdEdxReconUtils::GetPadGain warning!! indices out of range %d %d %d -- %d %d\n\n", det, icol, irow, calGainFactorROC->GetNcols(), calGainFactorROC->GetNrows() );  
175   }
176
177   return padgain;
178 }
179
180 Double_t AliTRDdEdxReconUtils::GetRNDClusterQ(AliTRDcluster *cl, const Double_t baseline)
181 {
182   //
183   //get cluter q from GetRawQ, apply baseline and Kr pad-calibration
184   //
185
186   const Int_t det     = cl->GetDetector();
187   const Int_t pad3col = cl->GetPadCol();
188   const Int_t padrow  = cl->GetPadRow();
189
190   Double_t rndqsum = 0;
191   for(Int_t ii=0; ii<7; ii++){
192     if(cl->GetSignals()[ii] < EPSILON){//bad pad marked by electronics
193       continue;
194     }
195
196     const Int_t icol = pad3col+(ii-3);
197     const Double_t padgain = GetPadGain(det, icol, padrow);
198     if(padgain<0){//indices out of range, pad3col near boundary case
199       continue;
200     }
201
202     const Double_t rndsignal = (cl->GetSignals()[ii] - baseline )/(AliTRDdEdxBaseUtils::IsPadGainOn()? padgain : 1);
203
204     //sum it anyway even if signal below baseline, as long as the total is positive
205     rndqsum += rndsignal;
206   }
207
208   return rndqsum;
209 }
210
211 Double_t AliTRDdEdxReconUtils::GetClusterQ(const Bool_t kinvq, const AliTRDseedV1 * seed, const Int_t itb)
212 {
213   //
214   //get cluster charge
215   //
216   Double_t dq = 0;
217   AliTRDcluster *cl = 0x0;
218       
219   const Double_t baseline = 10;
220
221   //GetRNDClusterQ(cl)>0 ensures that the total sum of q is above baseline*NsignalPhysical. 
222   cl = seed->GetClusters(itb);                    if(cl && GetRNDClusterQ(cl, baseline)>0 ) dq+= GetRNDClusterQ(cl, baseline);
223   cl = seed->GetClusters(itb+AliTRDseedV1::kNtb); if(cl && GetRNDClusterQ(cl, baseline)>0 ) dq+= GetRNDClusterQ(cl, baseline);
224
225   dq /= AliTRDdEdxBaseUtils::Getdldx(seed);
226   
227   dq /= AliTRDdEdxBaseUtils::QScale();
228       
229   if(kinvq){
230     if(dq>EPSILON){
231       dq = 1/dq;
232     }
233   }
234
235   return dq;
236 }
237
238 Int_t AliTRDdEdxReconUtils::GetArrayClusterQ(const Bool_t kinvq, TVectorD *arrayQ, TVectorD *arrayX, const AliTRDtrackV1 *trdtrack, Int_t timeBin0, Int_t timeBin1, Int_t tstep)
239 {
240   //
241   //return nclustter
242   //(if kinvq, return 1/q array), size of array must be larger than 31*6
243   //
244   if(!arrayQ || arrayQ->GetNrows()< (AliTRDseedV1::kNtb*AliTRDtrackV1::kNplane)){
245     printf("AliTRDdEdxReconUtils::GetArrayClusterQ error arrayQ null or size too small! %d\n", arrayQ? arrayQ->GetNrows() : -999); exit(1);
246   }
247   if(!arrayX || arrayX->GetNrows()< (AliTRDseedV1::kNtb*AliTRDtrackV1::kNplane)){
248     printf("AliTRDdEdxReconUtils::GetArrayClusterQ error arrayX null or size too small! %d\n", arrayX? arrayX->GetNrows() : -999); exit(1);
249   }
250
251   const Int_t mintb = 0;
252   const Int_t maxtb = AliTRDseedV1::kNtb-1;
253   if(timeBin0<mintb) timeBin0=mintb;
254   if(timeBin1>maxtb) timeBin1=maxtb;
255   if(tstep<=0) tstep=1;
256
257   //============
258   Int_t tbN=0;
259   Double_t tbQ[200];
260   Int_t tbBin[200];
261     
262   for(Int_t ichamber=0; ichamber < AliTRDtrackV1::kNplane; ichamber++){
263     const AliTRDseedV1 * seed = trdtrack->GetTracklet(ichamber);
264     if(!seed)
265       continue;
266     
267     const Int_t det = seed->GetDetector();
268
269     for(Int_t itb=timeBin0; itb<=timeBin1; itb+=tstep){
270       const Double_t dq = GetClusterQ(kinvq, seed, itb);
271       if(dq<EPSILON)
272         continue;
273
274       const Int_t gtb = det * AliTRDseedV1::kNtb + itb;
275
276       tbQ[tbN]=dq;
277       tbBin[tbN]=gtb;
278       tbN++;
279     }
280   }
281
282   Int_t ncls = 0;
283   for(Int_t iq=0; iq<tbN; iq++){
284     if(tbQ[iq]<EPSILON)
285       continue;
286
287     (*arrayQ)[ncls] = tbQ[iq];
288     (*arrayX)[ncls] = tbBin[iq];
289
290     ncls++;
291   }
292
293   static Int_t kprint = 100;
294   if(kprint<0){
295     printf("\nAliTRDdEdxReconUtils::GetArrayClusterQ raw cluster-Q\n");
296     for(Int_t iq=0; iq<ncls; iq++){
297       const Int_t ichamber =  AliTRDdEdxBaseUtils::ToLayer((*arrayX)[iq]);
298       const AliTRDseedV1 * seed = trdtrack->GetTracklet(ichamber);
299       if(!seed){
300         printf("error seed null!!\n"); exit(1);
301       }
302       const Double_t rawq =  (*arrayQ)[iq] * 45. * AliTRDdEdxBaseUtils::Getdldx(seed);
303       printf("esdid=%d; chamber=%d; timebin=%d; rawq= %.3f; myq[%d]= %e;\n", trdtrack->GetESDid(), ichamber, AliTRDdEdxBaseUtils::ToTimeBin((*arrayX)[iq]), rawq, iq, (*arrayQ)[iq]);
304     }
305     printf("\n");
306   }
307   kprint++;
308
309   return ncls;
310 }
311
312 Int_t AliTRDdEdxReconUtils::UpdateArrayX(const Int_t ncls, TVectorD* arrayX)
313 {
314   //
315   //arrayX det*Ntb+itb -> itb
316   //
317
318   TVectorD countChamber(6);
319   for(Int_t ii = 0; ii<ncls; ii++){
320     const Int_t xx = (Int_t)(*arrayX)[ii];
321     const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(xx);
322     
323     const Double_t ich = AliTRDgeometry::GetLayer(idet);
324     const Double_t itb = AliTRDdEdxBaseUtils::ToTimeBin(xx);
325     (*arrayX)[ii] = ich*AliTRDseedV1::kNtb+itb;
326
327     countChamber[ich] = 1;
328   }
329
330   const Double_t nch = countChamber.Sum();
331   return (Int_t) nch;
332 }
333
334