TRD/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-12
  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::GetAngularCorrection(const AliTRDseedV1 *seed)
 150 {
 151   //
 152   //return angular normalization factor
 153   //
 154
 155   return TMath::Sqrt(1+seed->GetYref(1)*seed->GetYref(1)+seed->GetZref(1)*seed->GetZref(1));
 156 }
 157
 158 Double_t AliTRDdEdxReconUtils::GetPadGain(const Int_t det, const Int_t icol, const Int_t irow)
 159 {
 160   //
 161   //get pad calibration
 162   //
 163   AliTRDcalibDB* calibration = AliTRDcalibDB::Instance();
 164   if (!calibration) {
 165     printf("AliTRDdEdxReconUtils::GetPadCalib No AliTRDcalibDB instance available\n"); exit(1);
 166   }
 167   AliTRDCalROC * calGainFactorROC = calibration->GetGainFactorROC(det);
 168   if(!calGainFactorROC){
 169     printf("AliTRDdEdxReconUtils::GetPadCalib no calGainFactorROC!\n"); exit(1);
 170   }
 171
 172   Double_t padgain = -999;
 173   if( icol >= 0 &&
 174       icol < calGainFactorROC->GetNcols() &&
 175       irow >=0 &&
 176       irow < calGainFactorROC->GetNrows()){
 177     padgain = calGainFactorROC->GetValue(icol, irow);
 178     if(padgain<EPSILON){
 179       printf("AliTRDdEdxReconUtils::GetPadGain padgain 0! %f %f -- %d %d %d -- %d %d\n", padgain, EPSILON, det, icol, irow, calGainFactorROC->GetNcols(), calGainFactorROC->GetNrows()); exit(1);
 180     }
 181   }
 182   else{
 183     //printf("\nAliTRDdEdxReconUtils::GetPadGain warning!! indices out of range %d %d %d -- %d %d\n\n", det, icol, irow, calGainFactorROC->GetNcols(), calGainFactorROC->GetNrows() );
 184   }
 185
 186   return padgain;
 187 }
 188
 189 Double_t AliTRDdEdxReconUtils::GetRNDClusterQ(AliTRDcluster *cl)
 190 {
 191   //
 192   //get cluter q from GetRawQ, apply baseline and Kr pad-calibration
 193   //
 194
 195   const Int_t det     = cl->GetDetector();
 196   const Int_t pad3col = cl->GetPadCol();
 197   const Int_t padrow  = cl->GetPadRow();
 198
 199   const Double_t baseline = 10;
 200
 201   Double_t rndqsum = 0;
 202   for(Int_t ii=0; ii<7; ii++){
 203     if(cl->GetSignals()[ii] < EPSILON){//bad pad marked by electronics
 204       continue;
 205     }
 206
 207     const Int_t icol = pad3col+(ii-3);
 208     const Double_t padgain = GetPadGain(det, icol, padrow);
 209     if(padgain<0){//indices out of range, pad3col near boundary case
 210       continue;
 211     }
 212
 213     const Double_t rndsignal = (cl->GetSignals()[ii] - baseline)/(AliTRDdEdxBaseUtils::IsPadGainOn()? padgain : 1);
 214
 215     //sum it anyway even if signal below baseline, as long as the total is positive
 216     rndqsum += rndsignal;
 217   }
 218
 219   return rndqsum;
 220 }
 221
 222 Double_t AliTRDdEdxReconUtils::GetClusterQ(const Bool_t kinvq, const AliTRDseedV1 * seed, const Int_t itb)
 223 {
 224   //
 225   //get cluster charge
 226   //
 227   Double_t dq = 0;
 228   AliTRDcluster *cl = 0x0;
 229
 230   //GetRNDClusterQ(cl)>0 ensures that the total sum of q is above baseline*NsignalPhysical.
 231   cl = seed->GetClusters(itb);                    if(cl && GetRNDClusterQ(cl)>0 ) dq+= GetRNDClusterQ(cl);//cl->GetRawQ();
 232   cl = seed->GetClusters(itb+AliTRDseedV1::kNtb); if(cl && GetRNDClusterQ(cl)>0 ) dq+= GetRNDClusterQ(cl);//cl->GetRawQ();
 233
 234   dq /= GetAngularCorrection(seed);
 235
 236   dq /= AliTRDdEdxBaseUtils::QScale();
 237
 238   if(kinvq){
 239     if(dq>EPSILON){
 240       dq = 1/dq;
 241     }
 242   }
 243
 244   return dq;
 245 }
 246
 247 Int_t AliTRDdEdxReconUtils::GetArrayClusterQ(const Bool_t kinvq, TVectorD *arrayQ, TVectorD *arrayX, const AliTRDtrackV1 *trdtrack, Int_t timeBin0, Int_t timeBin1, Int_t tstep)
 248 {
 249   //
 250   //return nclustter
 251   //(if kinvq, return 1/q array), size of array must be larger than 31*6
 252   //
 253   if(!arrayQ || arrayQ->GetNrows()< (AliTRDseedV1::kNtb*AliTRDtrackV1::kNplane)){
 254     printf("AliTRDdEdxReconUtils::GetArrayClusterQ error arrayQ null or size too small! %d\n", arrayQ? arrayQ->GetNrows() : -999); exit(1);
 255   }
 256   if(!arrayX || arrayX->GetNrows()< (AliTRDseedV1::kNtb*AliTRDtrackV1::kNplane)){
 257     printf("AliTRDdEdxReconUtils::GetArrayClusterQ error arrayX null or size too small! %d\n", arrayX? arrayX->GetNrows() : -999); exit(1);
 258   }
 259
 260   const Int_t mintb = 0;
 261   const Int_t maxtb = AliTRDseedV1::kNtb-1;
 262   if(timeBin0<mintb) timeBin0=mintb;
 263   if(timeBin1>maxtb) timeBin1=maxtb;
 264   if(tstep<=0) tstep=1;
 265
 266   //============
 267   Int_t tbN=0;
 268   Double_t tbQ[200];
 269   Int_t tbBin[200];
 270
 271   for(Int_t ichamber=0; ichamber < AliTRDtrackV1::kNplane; ichamber++){
 272     const AliTRDseedV1 * seed = trdtrack->GetTracklet(ichamber);
 273     if(!seed)
 274       continue;
 275
 276     const Int_t det = seed->GetDetector();
 277
 278     for(Int_t itb=timeBin0; itb<=timeBin1; itb+=tstep){
 279       const Double_t dq = GetClusterQ(kinvq, seed, itb);
 280       if(dq<EPSILON)
 281         continue;
 282
 283       const Int_t gtb = det * AliTRDseedV1::kNtb + itb;
 284
 285       tbQ[tbN]=dq;
 286       tbBin[tbN]=gtb;
 287       tbN++;
 288     }
 289   }
 290
 291   Int_t ncls = 0;
 292   for(Int_t iq=0; iq<tbN; iq++){
 293     if(tbQ[iq]<EPSILON)
 294       continue;
 295
 296     (*arrayQ)[ncls] = tbQ[iq];
 297     (*arrayX)[ncls] = tbBin[iq];
 298
 299     ncls++;
 300   }
 301
 302   static Int_t kprint = 100;
 303   if(kprint<0){
 304     printf("\nAliTRDdEdxReconUtils::GetArrayClusterQ raw cluster-Q\n");
 305     for(Int_t iq=0; iq<ncls; iq++){
 306       const Int_t ichamber =  AliTRDdEdxBaseUtils::ToLayer((*arrayX)[iq]);
 307       const AliTRDseedV1 * seed = trdtrack->GetTracklet(ichamber);
 308       if(!seed){
 309         printf("error seed null!!\n"); exit(1);
 310       }
 311       const Double_t rawq =  (*arrayQ)[iq] * 45. * GetAngularCorrection(seed);
 312       printf("esdid=%d; chamber=%d; timebin=%d; rawq= %.3f; myq[%d]= %e;\n", trdtrack->GetESDid(), ichamber, AliTRDdEdxBaseUtils::ToTimeBin((*arrayX)[iq]), rawq, iq, (*arrayQ)[iq]);
 313     }
 314     printf("\n");
 315   }
 316   kprint++;
 317
 318   return ncls;
 319 }
 320
 321 Int_t AliTRDdEdxReconUtils::UpdateArrayX(const Int_t ncls, TVectorD* arrayX)
 322 {
 323   //
 324   //arrayX det*Ntb+itb -> itb
 325   //
 326
 327   TVectorD countChamber(6);
 328   for(Int_t ii = 0; ii<ncls; ii++){
 329     const Int_t xx = (Int_t)(*arrayX)[ii];
 330     const Int_t idet = AliTRDdEdxBaseUtils::ToDetector(xx);
 331
 332     const Double_t ich = AliTRDgeometry::GetLayer(idet);
 333     const Double_t itb = AliTRDdEdxBaseUtils::ToTimeBin(xx);
 334     (*arrayX)[ii] = ich*AliTRDseedV1::kNtb+itb;
 335
 336     countChamber[ich] = 1;
 337   }
 338
 339   const Double_t nch = countChamber.Sum();
 340   return (Int_t) nch;
 341 }
 342
 343