Coverity fix
[u/mrichter/AliRoot.git] / TRD / TRDbase / AliTRDdEdxReconUtils.cxx
CommitLineData
6951a056 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
6fa94550 57#define EPSILON 1e-8
6951a056 58
59Int_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
91Double_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
123Double_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
6951a056 149Double_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
ca66da8f 180Double_t AliTRDdEdxReconUtils::GetRNDClusterQ(AliTRDcluster *cl, const Double_t baseline)
6951a056 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
6951a056 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
ca66da8f 202 const Double_t rndsignal = (cl->GetSignals()[ii] - baseline )/(AliTRDdEdxBaseUtils::IsPadGainOn()? padgain : 1);
6951a056 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
211Double_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
ca66da8f 219 const Double_t baseline = 10;
220
6951a056 221 //GetRNDClusterQ(cl)>0 ensures that the total sum of q is above baseline*NsignalPhysical.
ca66da8f 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);
6951a056 224
ca66da8f 225 dq /= AliTRDdEdxBaseUtils::Getdldx(seed);
6951a056 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
238Int_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 }
ca66da8f 302 const Double_t rawq = (*arrayQ)[iq] * 45. * AliTRDdEdxBaseUtils::Getdldx(seed);
6951a056 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
312Int_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