StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StdEdxY2Maker.cxx
1 // $Id: StdEdxY2Maker.cxx,v 1.101 2021/05/10 16:54:45 fisyak Exp $
2 //#define __NEGATIVE_ONLY__
3  #ifndef __NEGATIVE_ONLY__
4  #define __NEGATIVE_AND_POSITIVE__
5  #endif /* ! __NEGATIVE_ONLY__ */
6  #define __BEST_VERTEX__
7 #ifdef __TFG__VERSION__
8 //#define CompareWithToF
9 //#define __CHECK_LargedEdx__
10 //#define __TEST_DX__
11  #define __SpaceCharge__
12 //#define __LogProb__
13 //#define __DEBUG_dEdx__
14  #define __DEBUG_dNdx__
15 //#define __ADD_PROB__
16 //#define __BENCHMARKS__DOFIT_ZN__
17  #define __FIT_PULLS__
18  #define __CHECK_RDOMAP_AND_VOLTAGE__
19  #ifdef __CHECK_RDOMAP_AND_VOLTAGE__
20  #include "TProfile3D.h"
21  #endif /* __CHECK_RDOMAP_AND_VOLTAGE__ */
22 //#define __dZdY_dXdY__
23 //#define __Pad_Tmbk__
24 #endif /* __TFG__VERSION__ */
25 #include <Stiostream.h>
26 #include "StdEdxY2Maker.h"
27 #include "StTpcdEdxCorrection.h"
28 #include <vector>
29 // ROOT
30 #include "TMinuit.h"
31 #include "TSystem.h"
32 #include "TMath.h"
33 #include "TH2.h"
34 #include "TH3.h"
35 #include "THnSparse.h"
36 #include "TF1.h"
37 #include "TStyle.h"
38 #include "TProfile.h"
39 #include "TProfile2D.h"
40 #include "TTree.h"
41 #include "TChain.h"
42 #include "TFile.h"
43 #include "TNtuple.h"
44 #include "TCanvas.h"
45 #include "TRandom.h"
46 #include "TClonesArray.h"
47 #include "TArrayI.h"
48 #include "TArrayD.h"
49 #ifdef __BENCHMARKS__DOFIT_ZN__
50 #include "TBenchmark.h"
51 #endif /* __BENCHMARKS__DOFIT_ZN__ */
52 // StUtilities
53 #include "StMagF.h"
54 #include "StDbUtilities/StMagUtilities.h"
55 #include "StMessMgr.h"
56 #include "StBichsel/Bichsel.h"
57 #include "StBichsel/StdEdxModel.h"
58 #include "StBichsel/StdEdxPull.h"
59 #include "StDetectorId.h"
60 #include "StDedxMethod.h"
61 // StarClassLibrary
62 #include "SystemOfUnits.h"
63 #ifndef ST_NO_NAMESPACES
64 using namespace units;
65 #endif
66 #include "StarClassLibrary/StParticleTypes.hh"
67 
68 // StDb
69 #include "StDbUtilities/StTpcCoordinateTransform.hh"
70 #include "StDbUtilities/StCoordinates.hh"
71 #include "StTpcDb/StTpcDb.h"
72 // StEvent
73 #include "StEventTypes.h"
74 #include "StProbPidTraits.h"
75 #ifdef CompareWithToF
76 #include "StBTofPidTraits.h"
77 #endif
78 #include "StTpcDedxPidAlgorithm.h"
79 #include "StDetectorDbMaker/St_tss_tssparC.h"
80 #include "StDetectorDbMaker/St_tpcAnodeHVavgC.h"
81 #include "StDetectorDbMaker/St_TpcAvgCurrentC.h"
82 #include "StDetectorDbMaker/St_TpcAvgPowerSupplyC.h"
83 #include "StDetectorDbMaker/St_trigDetSumsC.h"
84 #include "StDetectorDbMaker/StDetectorDbTpcRDOMasks.h"
85 #include "StPidStatus.h"
86 #include "dEdxHist.h"
87 #if defined(__CHECK_LargedEdx__) || defined( __DEBUG_dEdx__) || defined( __DEBUG_dNdx__)
88 #include "tables/St_g2t_track_Table.h"
89 #endif
90 const static Int_t tZero= 19950101;
91 static Int_t tMin = 20000101;
92 static Int_t tMax = 20220101;
93 const static TDatime t0(tZero,0);
94 const static Int_t timeOffSet = t0.Convert();
95 static Double_t tpcTime = -1;
96 Int_t StdEdxY2Maker::NdEdx = 0;
97 dEdxY2_t *StdEdxY2Maker::CdEdx = 0;
98 dEdxY2_t *StdEdxY2Maker::FdEdx = 0;
99 dEdxY2_t *StdEdxY2Maker::dEdxS = 0;
100 static Int_t numberOfSectors = 0;
101 static Int_t numberOfTimeBins = 0;
102 static Int_t NumberOfChannels = 8;
103 TGraph *StdEdxY2Maker::fdNdxGraph[3] = {0};
104 
105 //const static Double_t pMomin = 0.35; // range for dE/dx calibration
106 //const static Double_t pMomax = 0.75;
107 const Double_t pMoMIP = 0.526; // MIP from Heed bg = 3.77 => p_pion = 0.526
108 const Double_t pMomin = pMoMIP - 0.05; // 0.45;
109 const Double_t pMomax = pMoMIP + 0.05; // 0.55;
110 
111 Double_t StdEdxY2Maker::bField = 0;
112 Bool_t StdEdxY2Maker::fUsedNdx = kFALSE;
113 TH2F *StdEdxY2Maker::fIntegratedAdc = 0;
114 //______________________________________________________________________________
115 // QA histograms
116 const static Int_t fNZOfBadHits = 10 + StTpcdEdxCorrection::kTpcAllCorrections;
117 static TH1F **fZOfBadHits = 0;
118 static TH1F *fZOfGoodHits = 0;
119 static TH1F *fPhiOfGoodHits = 0;
120 static TH1F *fPhiOfBadHits = 0;
121 static TH1F *fTracklengthInTpcTotal = 0;
122 static TH1F *fTracklengthInTpc = 0;
123 static TH2F *fPadTbkAll = 0;
124 static TH2F *fPadTbkBad = 0;
125 #ifdef __SpaceCharge__
126 static TH2F *AdcSC = 0, *AdcOnTrack = 0, *dEOnTrack = 0;
127 #endif
128 #ifdef __CHECK_RDOMAP_AND_VOLTAGE__
129 static TH3F *AlivePads = 0;
130 static TProfile3D *ActivePads = 0;
131 #endif
132 #ifdef __BEST_VERTEX__
133 static TH3F *PVxyz = 0, *PVxyzC = 0;
134 static TH2F *EtaVspT[2][2] = {0}; // Global and Primary, Positive and Negative
135 static TH2F *EtaVspTC[2] = {0}; // Positive and Negative global track used for calibration
136 #endif /* __BEST_VERTEX__ */
137 //______________________________________________________________________________
138 ClassImp(StdEdxY2Maker);
139 //_____________________________________________________________________________
140 StdEdxY2Maker::StdEdxY2Maker(const char *name): StMaker(name), m_Mask(-1) {
141  memset (beg, 0, end-beg);
142  SETBIT(m_Mode,kPadSelection);
143  m_Minuit = new TMinuit(2);
144  SetAttr("tMin",tMin);
145  SetAttr("tMax",tMax);
146 }
147 //_____________________________________________________________________________
148 Int_t StdEdxY2Maker::Init(){
149  if (IAttr("EmbeddingShortCut")) {
150  m_Mode = 0;
151  SETBIT(m_Mode,kEmbedding);
152  SETBIT(m_Mode,kPadSelection);
153  SETBIT(m_Mask,StTpcdEdxCorrection::kTpcLast);
154  }
155  if (Debug()) {
156  if (! TESTBIT(m_Mode, kOldClusterFinder)) {
157  LOG_WARN << "StdEdxY2Maker::Init use new Cluster Finder parameterization" << endm;
158  } else {
159  LOG_WARN << "StdEdxY2Maker::Init use old Cluster Finder parameterization" << endm;
160  }
161  if (TESTBIT(m_Mode, kPadSelection)) {
162  LOG_WARN << "StdEdxY2Maker::Init Pad Selection is ON" << endm;
163  }
164  if (TESTBIT(m_Mode, kDoNotCorrectdEdx)) {
165  LOG_WARN << "StdEdxY2Maker::Init Don't Correct dEdx" << endm;
166  }
167  if (TESTBIT(m_Mode, kEmbedding)) {
168  LOG_WARN << "StdEdxY2Maker::Init This is embedding run" << endm;
169  }
170  }
171  gMessMgr->SetLimit("StdEdxY2Maker:: mismatched Sector",20);
172  gMessMgr->SetLimit("StdEdxY2Maker:: pad/TimeBucket out of range:",20);
173  gMessMgr->SetLimit("StdEdxY2Maker:: Helix Prediction",20);
174  gMessMgr->SetLimit("StdEdxY2Maker:: Coordinates",20);
175  gMessMgr->SetLimit("StdEdxY2Maker:: Prediction",20);
176  gMessMgr->SetLimit("StdEdxY2Maker:: NdEdx",20);
177  gMessMgr->SetLimit("StTpcdEdxCorrection:: Illegal time for scalers",20);
178  return StMaker::Init();
179 }
180 //_____________________________________________________________________________
181 Int_t StdEdxY2Maker::InitRun(Int_t RunNumber){
182  tMin = IAttr("tMin");
183  tMax = IAttr("tMax");
184  static Int_t DoOnce = 0;
185  if (!gStTpcDb) {
186  cout << "Database Missing! Can't initialize StdEdxY2Maker" << endl;
187  return kStFatal;
188  }
189  // TPG parameters
190  numberOfSectors = gStTpcDb->Dimensions()->numberOfSectors();
191  numberOfTimeBins = gStTpcDb->Electronics()->numberOfTimeBins();
192  SafeDelete(m_TpcdEdxCorrection);
193  m_TpcdEdxCorrection = new StTpcdEdxCorrection(m_Mask, Debug());
194 
195  if (! DoOnce) {
196  DoOnce = 1;
197  if (! IAttr("SkipdNdx")) {
198  if ((GetDate() > 20171201 && m_TpcdEdxCorrection->IsFixedTarget()) ||
199  (GetDate() > 20181201)) fUsedNdx = kTRUE; // use dN/dx for fixed target for Run XVIII and year >= XIX
200  }
201  if (TESTBIT(m_Mode, kCalibration)) {// calibration mode
202  if (Debug()) LOG_WARN << "StdEdxY2Maker::InitRun Calibration Mode is On (make calibration histograms)" << endm;
203  TFile *f = GetTFile();
204  if (f) {
205  f->cd();
206  if ((TESTBIT(m_Mode, kGASHISTOGRAMS))) {
207  if (Debug()) LOG_WARN << "StdEdxY2Maker::InitRun Gas Histograms is ON" << endm;
208  TrigHistos();
209  }
210  Histogramming();
211  if (TESTBIT(m_Mode, kXYZcheck)) XyzCheck();
212  }
213  }
214  QAPlots(0);
215  // Switch between usage prediction dependence of log2(dX), new option is not to use it in predicetion due to PicoDst
216  fUsedx2 = kFALSE;
217  if (m_TpcdEdxCorrection->Correction(StTpcdEdxCorrection::kTpcLengthCorrectionMD2) ||
218  m_TpcdEdxCorrection->Correction(StTpcdEdxCorrection::kTpcLengthCorrectionMDF) ||
219  m_TpcdEdxCorrection->Correction(StTpcdEdxCorrection::kTpcLengthCorrection )) {
220  fUsedx2 = kTRUE;
221  }
222  if (m_TpcdEdxCorrection->Correction(StTpcdEdxCorrection::kTpcLengthCorrectionMDN)) {
223  fUsedx2 = kFALSE;
224  }
225  LOG_WARN << "StdEdxY2Maker::InitRun Force ";
226  if (fUsedx2) {
227  LOG_WARN << "to USE";
228  } else {
229  LOG_WARN << "NOT to USE";
230  }
231  LOG_WARN << " dx2L in dE/dx predictions "<< endm;
232  }
233  St_tpcPadGainT0C::instance(); // activate extra gain corrections for tpx
234  St_itpcPadGainT0C::instance(); // activate extra gain corrections for iTPC
235  return kStOK;
236 }
237 //_______________________________________________________________________________
238 Double_t StdEdxY2Maker::gaus2(Double_t *x, Double_t *p) {
239  Double_t NormL = p[0];
240  Double_t mu = p[1];
241  Double_t muP = mu + p[4];
242  Double_t sigma = p[2];
243  Double_t sigmaP = TMath::Sqrt(sigma*sigma + 0.101741*0.101741);
244  Double_t phi = p[3];
245  Double_t frac = TMath::Sin(phi);
246  frac *= frac;
247  return TMath::Exp(NormL)*((1 - frac)*TMath::Gaus(x[0],mu ,sigma ,kTRUE) +
248  frac *TMath::Gaus(x[0],muP,sigmaP,kTRUE));
249 }
250 //________________________________________________________________________________
251 TF1 *StdEdxY2Maker::Gaus2() {
252  TF1 *f = new TF1("Gaus2",gaus2,-3,3,5);
253  f->SetParName(0,"NormL"); f->SetParLimits(0,-10.,10.);
254  f->SetParName(1,"mu"); f->SetParLimits(1,-0.5,0.5);
255  f->SetParName(2,"sigma"); f->SetParLimits(2, 0.2,0.5);
256  f->SetParName(3,"phiP"); f->SetParLimits(3, 0.0,TMath::Pi()/4);
257  f->SetParName(4,"muP");
258  f->SetParameters(10,0,0.3,0.1,1.315);
259  // f->FixParameter(4,1.425);
260  return f;
261 }
262 //_____________________________________________________________________________
263 Int_t StdEdxY2Maker::Finish() {
264  // static Double_t slope = 1.7502e-6;// slope from Blair 1/( O2 in ppm., cm )
265  SafeDelete(m_Minuit);
266  if (m_TpcdEdxCorrection && m_TpcdEdxCorrection->TestBit(kCanDelete)) delete m_TpcdEdxCorrection;
267  m_TpcdEdxCorrection = 0;
268  return StMaker::Finish();
269 }
270 //_____________________________________________________________________________
271 void StdEdxY2Maker::AddEdxTraits(StTrack *tracks[2], dst_dedx_st &dedx){
272  for (Int_t l = 0; l < 2; l++) {
273  if (tracks[l]) {
274  StDedxPidTraits *trait = new StDedxPidTraits(dedx);
275  tracks[l]->addPidTraits(trait);
276  }
277  }
278 }
279 //_____________________________________________________________________________
280 Int_t StdEdxY2Maker::Make(){
281 #ifdef __BENCHMARKS__DOFIT_ZN__
282  TBenchmark myBenchmark;
283 #endif /* __BENCHMARKS__DOFIT_ZN__ */
284  static Bool_t ForcedX = IAttr("ForcedX");
285  tpcTime = GetDateTime().Convert() - timeOffSet;
286  static StTpcLocalSectorCoordinate localSect[4];
287  static StTpcPadCoordinate PadOfTrack, Pad;
288  static StTpcLocalSectorDirection localDirectionOfTrack;
289  static StThreeVectorD xyz[4];
290  static StThreeVectorD dirG;
291  static Double_t s[2], s_in[2], s_out[2], w[2], w_in[2], w_out[2], dx, AdcI, dxC;
292 #ifdef __dZdY_dXdY__
293  static Double dZdY, dXdY;
294 #endif
295  enum {kNdEdxMax = 300};
296  static dEdxY2_t CdEdxT[3*kNdEdxMax];//,FdEdxT[kNdEdxMax],dEdxST[kNdEdxMax];
297  static Int_t sectorMin = 1, sectorMax = 24;
298  CdEdx = CdEdxT;
299  FdEdx = CdEdxT + kNdEdxMax; // sorted by row
300  dEdxS = CdEdxT + 2*kNdEdxMax; // sorted by dE/dx
301  St_tpcGas *tpcGas = m_TpcdEdxCorrection->tpcGas();
302  if (TESTBIT(m_Mode, kCalibration) && tpcGas) TrigHistos(1);
303  StTpcCoordinateTransform transform(gStTpcDb);
304  StEvent* pEvent = dynamic_cast<StEvent*> (GetInputDS("StEvent"));
305  if (!pEvent) {
306  LOG_INFO << "StdEdxY2Maker: no StEvent " << endm;
307  return kStOK; // if no event, we're done
308  }
309  if (pEvent->runInfo()) bField = pEvent->runInfo()->magneticField()*kilogauss;
310  if (TMath::Abs(bField) < 1.e-5*kilogauss) return kStOK;
311  UInt_t NoPV = pEvent->numberOfPrimaryVertices();
312  if (! NoPV) return kStOK;
313 #ifdef __BEST_VERTEX__
314  const StBTofCollection* tof = pEvent->btofCollection();
315  StPrimaryVertex *pVbest = pEvent->primaryVertex(0);
316  Double_t VpdZ = -300;
317  if (tof) {
318  if (tof->tofHeader()) VpdZ = tof->tofHeader()->vpdVz();
319  }
320  if (m_TpcdEdxCorrection->IsFixedTarget()) VpdZ = 200;
321  Double_t dZbest = 999;
322  for (UInt_t ipr = 0; ipr < NoPV; ipr++) {
323  StPrimaryVertex *pVertex = pEvent->primaryVertex(ipr);
324  if (! pVertex) continue;
325  if (PVxyz) PVxyz->Fill( pVertex->position().x(), pVertex->position().y(), pVertex->position().z());
326  if (TMath::Abs(VpdZ) < 250) {
327  Double_t zTPC = pVertex->position().z();
328  Double_t dZ = TMath::Abs(zTPC-VpdZ);
329  if (dZ < dZbest) {
330  dZbest = dZ;
331  pVbest = pVertex;
332  }
333  }
334  }
335  if (dZbest < 999 && dZbest > 3.0) pVbest = 0;
336  if (pVbest &&PVxyzC) PVxyzC->Fill( pVbest->position().x(), pVbest->position().y(), pVbest->position().z());
337 #endif /* __BEST_VERTEX__ */
338  // no of tpc hits
339  Int_t TotalNoOfTpcHits = 0;
340  Int_t NoOfTpcHitsUsed = 0;
341  StTpcHitCollection* TpcHitCollection = pEvent->tpcHitCollection();
342  if (! TpcHitCollection) {
343  LOG_INFO << "StdEdxY2Maker: no TpcHitCollection " << endm;
344  return kStOK; // if no event, we're done
345  }
346  TotalNoOfTpcHits = TpcHitCollection->numberOfHits();
347  StSPtrVecTrackNode& trackNode = pEvent->trackNodes();
348  UInt_t nTracks = trackNode.size();
349  for (UInt_t i=0; i < nTracks; i++) {
350  StTrackNode *node = trackNode[i];
351  if (!node) continue;
352  StGlobalTrack *gTrack = dynamic_cast<StGlobalTrack *>(node->track(global));
353  if (gTrack && gTrack->bad()) {gTrack = 0;}
354  if (! gTrack || gTrack->flag() <= 0) continue;
355  StPrimaryTrack *pTrack = 0;
356  if (gTrack) {
357  pTrack = static_cast<StPrimaryTrack*>(node->track(primary));
358  }
359  if (pTrack && pTrack->bad()) {pTrack = 0;}
360  StTrack *track = 0;
361  StTrack *tracks[2] = {gTrack, pTrack};
362  Int_t sCharge = 0; // positive
363  if (gTrack->geometry()->charge() < 0) sCharge = 1;// negative
364  Int_t qB = (sCharge+1)%2;
365  if (bField > 0) qB = (qB + 1)%2; // swap for Full Field
366 #ifdef __BEST_VERTEX__
367  if (TESTBIT(m_Mode, kCalibration)) {// calibration mode
368  for (Int_t l = 0; l < 2; l++) {
369  track = tracks[l];
370  if (track) {
371  StThreeVectorD g3 = track->geometry()->momentum(); // p of global track
372  EtaVspT[l][sCharge]->Fill(TMath::Log10(g3.perp()), g3.pseudoRapidity());
373  }
374  }
375  }
376 #endif /* __BEST_VERTEX__ */
377  StPhysicalHelixD helixO = gTrack->outerGeometry()->helix();
378  StPhysicalHelixD helixI = gTrack->geometry()->helix();
379  StThreeVectorD g3 = gTrack->geometry()->momentum(); // p of global track
380  Double_t etaG = g3.pseudoRapidity();
381  if (Debug() > 1) {
382  cout << "Track:" << i
383  << "\ttype " << gTrack->type()
384  << "\tvertex " << gTrack->vertex()
385  << "\tkey " << gTrack->key()
386  << "\tflag " << gTrack->flag()
387  << "\tencodedMethod " << gTrack->encodedMethod()
388  << "\timpactParameter " << gTrack->impactParameter()
389  << "\tlength " << gTrack->length()
390  << "\tEtaG " << etaG
391  << "\tnumberOfPossiblePoints " << gTrack->numberOfPossiblePoints() << endl;
392  cout << "pxyzI: " << gTrack->geometry()->momentum() << "\tmag " << gTrack->geometry()->momentum().mag() << endl;
393  cout << "pxyzO: " << gTrack->outerGeometry()->momentum() << "\tmag " << gTrack->outerGeometry()->momentum().mag() << endl;
394  cout << "start Point: " << helixI.at(0) << endl;
395  cout << "end Point: " << helixO.at(0) << endl;
396  }
397  // StPtrVecHit hvec = gTrack->detectorInfo()->hits(kTpcId);
398  StPtrVecHit hvec = gTrack->detectorInfo()->hits();
399  // if no hits than make only histograms. Works if kDoNotCorrectdEdx mode is set
400  if (! TESTBIT(m_Mode, kDoNotCorrectdEdx)) {
401  // clean old PiD traits
402  for (Int_t l = 0; l < 2; l++) {
403  track = tracks[l];
404  if (track) {
405  StSPtrVecTrackPidTraits &traits = track->pidTraits();
406  UInt_t size = traits.size();
407  if (size) {
408  for (UInt_t i = 0; i < size; i++) {
409  StDedxPidTraits *pid = dynamic_cast<StDedxPidTraits*>(traits[i]);
410  if (! pid) continue;
411  if (pid->detector() != kTpcId) continue;
412  traits[i]->makeZombie(1);
413  }
414  }
415  }
416  }
417  }
418  if (hvec.size() && ! TESTBIT(m_Mode, kDoNotCorrectdEdx)) {
419  Int_t Id = gTrack->key();
420  Int_t NoFitPoints = gTrack->fitTraits().numberOfFitPoints(kTpcId);
421  NdEdx = 0;
422  Double_t TrackLength70 = 0;
423  Double_t TrackLength = 0;
424  Double_t TrackLengthTotal = 0;
425  for (UInt_t j=0; j<hvec.size(); j++) {// hit loop
426  if (hvec[j]->detector() != kTpcId) continue;
427  StTpcHit *tpcHit = static_cast<StTpcHit *> (hvec[j]);
428  if (! tpcHit) continue;
429  if (Debug() > 1) {tpcHit->Print();}
430  if (! tpcHit->usedInFit()) {
431  BadHit(0,tpcHit->position());
432  continue;
433  } if ( tpcHit->flag()) {
434  BadHit(1,tpcHit->position());
435  continue;
436  }
437  Int_t sector = tpcHit->sector();
438  if (sector < sectorMin || sector > sectorMax) continue;
439  Int_t row = tpcHit->padrow();
440  Int_t pad = tpcHit->pad();
441  Int_t iRdo = StDetectorDbTpcRDOMasks::instance()->rdoForPadrow(sector,row,pad);
442  if ( ! StDetectorDbTpcRDOMasks::instance()->isOn(sector,iRdo)) continue;
443  if (! St_tpcAnodeHVavgC::instance()->livePadrow(sector,row)) continue;
444  xyz[3] = StThreeVectorD(tpcHit->position().x(),tpcHit->position().y(),tpcHit->position().z());
445  //________________________________________________________________________________
446  Float_t dX_TrackFit = tpcHit->dX();
447  Float_t dX_Helix = 0;
448  dx = dX_TrackFit;
449  AdcI = 0;
450 #ifdef __dZdY_dXdY__
451  dZdY = dXdY = 0;
452 #endif
453  static StGlobalDirection globalDirectionOfTrack;
454  Int_t iokCheck = 0;
455 #ifdef __TEST_DX__
456  static Bool_t TestdX = kTRUE;
457 #else
458  static Bool_t TestdX = kFALSE;
459 #endif /* __TEST_DX__ */
460  // use cluster position for precalculated dx
461  transform(xyz[3],localSect[3],sector,row);
462  transform(localSect[3],Pad);
463  while ((ForcedX || dX_TrackFit <= 0.0 || TestdX)) {
464  dX_Helix = -13;
465  StThreeVectorD middle = xyz[3];
466  StThreeVectorD upper(tpcHit->positionU().x(),tpcHit->positionU().y(),tpcHit->positionU().z());
467  StThreeVectorD lower(tpcHit->positionL().x(),tpcHit->positionL().y(),tpcHit->positionL().z());
468  StThreeVectorD dif = upper - lower;
469  StThreeVectorD normal = dif.unit();
470  StGlobalCoordinate globalOfTrack;
471  Double_t pad;
472 #if 0
473  StThreeVectorD &V = *&normal;
474  Double_t zd = sector <=12 ? 1: -1;
475  StThreeVectorD W = StThreeVectorD(0,0,zd);
476  StThreeVectorD U = V.cross(W);
477  StThreeVectorD D = dif.unit();
478  Double_t dY = 0;
479 #endif
480  // check that helix prediction is consistent with measurement
481  if (Propagate(middle,normal,helixI,helixO,xyz[0],dirG,s,w)) break;
482  if (Debug() > 1) {
483  cout << " Prediction:\t" << xyz[0]
484  << "\tat s=\t" << s[0] << "/" << s[1]
485  << "\tw = " << w[0] << "/" << w[1] << endl;
486  }
487  dif = xyz[3] - xyz[0];
488  if (dif.perp() > 2.0) {
489  if (Debug() > 1) {cout << "Prediction is to far from hit:\t" << xyz[3] << endl;}
490  break;
491  }
492  if (Propagate(upper,normal,helixI,helixO,xyz[1],dirG,s_out,w_out)) break;
493  if (Propagate(lower,normal,helixI,helixO,xyz[2],dirG,s_in ,w_in )) break;
494  dX_Helix = ((s_out[0] - s_in[0])*w[1] + (s_out[1] - s_in[1])*w[0]);
495  dif = xyz[1] - xyz[2];
496  // Check for Membernane
497  if (xyz[1].z() * xyz[2].z() < 0) {
498  Double_t dZ = TMath::Abs(xyz[1].z()) + TMath::Abs(xyz[2].z());
499  Double_t scaledX = 1;
500  if (xyz[1].z() * xyz[3].z() > 0) {
501  scaledX = TMath::Abs(xyz[1].z())/dZ;
502  } else if (xyz[2].z() * xyz[3].z() > 0) {
503  scaledX = TMath::Abs(xyz[2].z())/dZ;
504  }
505  static Int_t ibreak = 0;
506  if (Debug() > 1) {
507  cout << "Cross Membrane : upper " << xyz[1] << endl;
508  cout << " hit " << xyz[3] << endl;
509  cout << " lower " << xyz[2] << "\tscale dX = " << scaledX << endl;
510  }
511  dX_Helix *= scaledX;
512  ibreak++;
513  }
514  if (dX_Helix <= 0.0) {
515  if (Debug() > 1) {cout << "negative dX_Helix " << dX_Helix << endl;}
516  break;
517  }
518  // Consistency check
519  globalDirectionOfTrack = StGlobalDirection(dirG);
520  for (Int_t l = 0; l < 4; l++) {
521  globalOfTrack = StGlobalCoordinate(xyz[l].x(),xyz[l].y(),xyz[l].z());
522  transform(globalOfTrack,localSect[l],sector,row);
523  }
524  if (ForcedX || dX_TrackFit <= 0.0 )
525  tpcHit->setdX(dx);
526  transform(localSect[0],PadOfTrack);
527  transform(globalDirectionOfTrack,localDirectionOfTrack,sector,row);
528  transform(localSect[3],Pad);
529  if (sector != Pad.sector() || // ? && TMath::Abs(xyz[0].x()) > 20.0 ||
530  row != Pad.row()) {
531  LOG_WARN << "StdEdxY2Maker:: mismatched Sector "
532  << Pad.sector() << " / " << sector
533  << " Row " << Pad.row() << " / " << row
534  << "pad " << Pad.pad() << " TimeBucket :" << Pad.timeBucket()
535  << endm;
536  iokCheck++;
537  }
538  pad = tpcHit->pad();
539  if (pad == 0) pad = Pad.pad();
540  if (Pad.timeBucket() < 0 ||
541  Pad.timeBucket() >= numberOfTimeBins) {
542  LOG_WARN << "StdEdxY2Maker:: TimeBucket out of range: "
543  << Pad.timeBucket() << endm;
544  iokCheck++;
545  }
546  if (sector != PadOfTrack.sector() ||
547  row != PadOfTrack.row() ||
548  TMath::Abs(Pad.pad()-PadOfTrack.pad()) > 5) {
549  if (Debug() > 1) {
550  LOG_WARN << "StdEdxY2Maker:: Helix Prediction "
551  << "Sector = "
552  << PadOfTrack.sector() << "/"
553  << sector
554  << " Row = " << PadOfTrack.row() << "/"
555  << row
556  << " Pad = " << PadOfTrack.pad() << "/"
557  << Pad.pad()
558  << " from Helix is not matched with point/" << endm;;
559  LOG_WARN << "StdEdxY2Maker:: Coordinates Preiction: "
560  << xyz[0] << "/Hit " << tpcHit->position()
561  << endm;
562  }
563  iokCheck++;
564  }
565  if (iokCheck) {
566  dX_Helix = -13;
567  break;
568  }
569 #ifdef __dZdY_dXdY__
570  dY = D.dot(V);
571  if (TMath::Abs(dY) > 1e-7) {
572  dZdY = D.dot(W)/dY;
573  dXdY = D.dot(U)/dY;
574  } else {
575  dZdY = dXdY = 0;
576  }
577 #endif
578  if (Debug() > 1) {
579  cout << "Helix Prediction with dX = " << dX_Helix << endl;
580  }
581  break;
582  } // end of dx calculation
583  if (dX_Helix < 0) {
584  if (Debug() > 1) {
585  cout << "Helix Prediction Failed" << endl;
586  }
587  }
588  dx = tpcHit->dX();
589  if ((ForcedX || dX_TrackFit <= 0.0)) {
590  if (dX_Helix <= 0.0) continue;
591  if (Debug()) {
592  }
593  dx = dX_Helix;
594  if (ForcedX) tpcHit->setdX(dx);
595  } else {
596  dx = dX_TrackFit;
597  }
598  TrackLengthTotal += dx;
599  //________________________________________________________________________________
600  if (tpcHit->adc() <= 0) {
601  LOG_WARN << "StdEdxY2Maker:: adc : " << tpcHit->adc()
602  << " <= 0" << endm;
603  iokCheck++;
604  }
605  // if ((TESTBIT(m_Mode, kXYZcheck)) && (TESTBIT(m_Mode, kCalibration))) XyzCheck(&global, iokCheck);
606  if ((TESTBIT(m_Mode, kPadSelection)) && iokCheck) {BadHit(3, tpcHit->position()); continue;}
607  if ((TESTBIT(m_Mode, kPadSelection)) && (dx < 0.5 || dx > 25.)) {BadHit(4, tpcHit->position()); continue;}
608  // Corrections
609  CdEdx[NdEdx].Reset();
610  CdEdx[NdEdx].resXYZ[0] = localSect[3].position().x() - localSect[0].position().x();
611  CdEdx[NdEdx].resXYZ[1] = localSect[3].position().y() - localSect[0].position().y();
612  CdEdx[NdEdx].resXYZ[2] = localSect[3].position().z() - localSect[0].position().z();
613  CdEdx[NdEdx].DeltaZ = 5.2;
614  CdEdx[NdEdx].QRatio = -2;
615  CdEdx[NdEdx].QRatioA = -2.;
616  CdEdx[NdEdx].QSumA = 0;
617  CdEdx[NdEdx].sector = sector;
618  CdEdx[NdEdx].row = row;
619  CdEdx[NdEdx].pad = Pad.pad();
620  CdEdx[NdEdx].edge = CdEdx[NdEdx].pad;
621  Float_t NoPadsInRow = St_tpcPadConfigC::instance()->numberOfPadsAtRow(sector,row);
622  if (CdEdx[NdEdx].edge > 0.5*NoPadsInRow) CdEdx[NdEdx].edge -= 1 + NoPadsInRow;
623  CdEdx[NdEdx].xpad = 2*(CdEdx[NdEdx].pad - 0.5)/NoPadsInRow - 1.0;
624  CdEdx[NdEdx].xpadR = 2*(tpcHit->pad() - 0.5)/NoPadsInRow - 1.0;
625  CdEdx[NdEdx].qB = qB;
626  CdEdx[NdEdx].yrow = sector + 0.5*((row <= St_tpcPadConfigC::instance()->innerPadRows(sector)) ?
627  (row - St_tpcPadConfigC::instance()->innerPadRows(sector) - 0.5)/St_tpcPadConfigC::instance()->innerPadRows(sector) :
628  (row - St_tpcPadConfigC::instance()->innerPadRows(sector) - 0.5)/(St_tpcPadConfigC::instance()->numberOfRows(sector) - St_tpcPadConfigC::instance()->innerPadRows(sector)));
629  CdEdx[NdEdx].dX_TrackFit = dX_TrackFit;
630  CdEdx[NdEdx].dX_Helix = dX_Helix;
631  CdEdx[NdEdx].Npads = tpcHit->padsInHit();
632  CdEdx[NdEdx].Ntbks = tpcHit->timeBucketsInHit();
633  CdEdx[NdEdx].dCharge = 0;
634  CdEdx[NdEdx].rCharge= 0.5*m_TpcdEdxCorrection->Adc2GeV()*TMath::Pi()/4.*CdEdx[NdEdx].Npads*CdEdx[NdEdx].Ntbks;
635  if (TESTBIT(m_Mode, kEmbeddingShortCut) &&
636  (tpcHit->idTruth() && tpcHit->qaTruth() > 95)) CdEdx[NdEdx].lSimulated = tpcHit->idTruth();
637 #ifdef __dZdY_dXdY__
638  CdEdx[NdEdx].dZdY = dZdY;
639  CdEdx[NdEdx].dXdY = dXdY;
640 #endif
641  CdEdx[NdEdx].AdcI = AdcI;
642  dxC = dx;
643 #if 1
644  // Scale dX to full pad length
645  if (St_tpcPadConfigC::instance()->iTpc(sector)) {
646  Int_t io = 1;
647  if (row > St_tpcPadConfigC::instance()->innerPadRows(sector)) io = 0;
648  Double_t padlength = (io == 1) ?
649  St_tpcPadConfigC::instance()->innerSectorPadLength(sector) :
650  St_tpcPadConfigC::instance()->outerSectorPadLength(sector);
651  Double_t rowPitch = (io == 1) ?
652  St_tpcPadConfigC::instance()->innerSectorRowPitch1(sector) :
653  St_tpcPadConfigC::instance()->outerSectorRowPitch(sector);
654  dxC *= rowPitch/padlength;
655  }
656 #endif
657  CdEdx[NdEdx].dxC = dxC;
658  CdEdx[NdEdx].F.dE = tpcHit->charge();
659  CdEdx[NdEdx].F.dx = dxC;
660  CdEdx[NdEdx].xyz[0] = localSect[3].position().x();
661  CdEdx[NdEdx].xyz[1] = localSect[3].position().y();
662  CdEdx[NdEdx].xyz[2] = localSect[3].position().z();
663  Double_t maxPad = NoPadsInRow/2;
664  Double_t pitch = (row <= St_tpcPadConfigC::instance()->innerPadRows(sector)) ?
665  St_tpcPadConfigC::instance()->innerSectorPadPitch(sector) :
666  St_tpcPadConfigC::instance()->outerSectorPadPitch(sector);
667  Double_t PhiMax = TMath::ATan2(maxPad*pitch, St_tpcPadConfigC::instance()->radialDistanceAtRow(sector,row));
668  CdEdx[NdEdx].PhiR = TMath::ATan2(CdEdx[NdEdx].xyz[0],CdEdx[NdEdx].xyz[1])/PhiMax;
669  CdEdx[NdEdx].xyzD[0] = localDirectionOfTrack.position().x();
670  CdEdx[NdEdx].xyzD[1] = localDirectionOfTrack.position().y();
671  CdEdx[NdEdx].xyzD[2] = localDirectionOfTrack.position().z();
672  CdEdx[NdEdx].ZdriftDistance = localSect[3].position().z();
673  CdEdx[NdEdx].zG = tpcHit->position().z();
674  CdEdx[NdEdx].etaG = etaG;
675  Double_t pT2 = CdEdx[NdEdx].xyzD[0]*CdEdx[NdEdx].xyzD[0]+CdEdx[NdEdx].xyzD[1]*CdEdx[NdEdx].xyzD[1];
676  if (pT2 > 1e-14)
677  CdEdx[NdEdx].TanL = -CdEdx[NdEdx].xyzD[2]/TMath::Sqrt(pT2);
678  CdEdx[NdEdx].tpcTime = tpcTime;
679  if (St_trigDetSumsC::instance()) CdEdx[NdEdx].Zdc = St_trigDetSumsC::instance()->zdcX();
680  CdEdx[NdEdx].adc = tpcHit->adc();
681  Bool_t doIT = kTRUE;
682  if (TESTBIT(m_Mode,kEmbedding)) doIT = kFALSE;
683  if (fPadTbkAll) fPadTbkAll->Fill(CdEdx[NdEdx].Ntbks, CdEdx[NdEdx].Npads);
684  Int_t iok = m_TpcdEdxCorrection->dEdxCorrection(CdEdx[NdEdx],doIT);
685  if (iok) {
686  BadHit(10+iok, tpcHit->position());
687  if (fPadTbkBad) fPadTbkBad->Fill(CdEdx[NdEdx].Ntbks, CdEdx[NdEdx].Npads);
688  continue;
689  }
690  if (fZOfGoodHits) fZOfGoodHits->Fill(tpcHit->position().z());
691  if (fPhiOfGoodHits!= 0) fPhiOfGoodHits->Fill(TMath::ATan2(tpcHit->position().y(),tpcHit->position().x()));
692  if (NdEdx < kNdEdxMax) {
693  tpcHit->setCharge(CdEdx[NdEdx].F.dE);
694  TrackLength += CdEdx[NdEdx].F.dx;
695  NdEdx++;
696  NoOfTpcHitsUsed++;
697  }
698  if (NdEdx > NoFitPoints)
699  LOG_ERROR << "StdEdxY2Maker:: NdEdx = " << NdEdx
700  << "> NoFitPoints ="<< NoFitPoints << endm;
701  }
702  if (fTracklengthInTpcTotal) fTracklengthInTpcTotal->Fill(TrackLengthTotal);
703  if (fTracklengthInTpc) fTracklengthInTpc->Fill(TrackLength);
704  SortdEdx();
705  if (Debug() > 1) PrintdEdx(2);
706  Double_t I70 = 0, D70 = 0;
707  Double_t dXavLog2 = 1;
708  Double_t SumdEdX = 0;
709  Double_t SumdX = 0;
710  Int_t N70 = NdEdx - (int) (0.3*NdEdx + 0.5);
711  if (N70 > 1) {
712  Int_t k;
713  for (k = 0; k < N70; k++) {
714  I70 += dEdxS[k].F.dEdx;
715  D70 += dEdxS[k].F.dEdx*dEdxS[k].F.dEdx;
716  TrackLength70 += dEdxS[k].F.dx;
717  if (dEdxS[k].F.dx > 0) {
718  SumdEdX += dEdxS[k].F.dEdx;
719  SumdX += dEdxS[k].F.dEdx*TMath::Log2(dEdxS[k].F.dx);
720  }
721  }
722  I70 /= N70; D70 /= N70;
723  D70 = TMath::Sqrt(TMath::Abs(D70 - I70*I70));
724  D70 /= I70;
725 #ifdef __CHECK_LargedEdx__
726  static Double_t dEdxMin = 6e-6; // 6keV
727  static Int_t iBreak = 0;
728  if (I70 > dEdxMin) {
729  iBreak++;
730  gTrack->Print();
731  if (gTrack->idTruth() > 0) {
732  St_g2t_track *g2t_track = (St_g2t_track *) GetDataSet("geant/g2t_track");
733  if (g2t_track) g2t_track->Print(gTrack->idTruth()-1,1);
734  }
735  PrintdEdx(2);
736  }
737 #endif
738  if (SumdEdX > 0) dXavLog2 = SumdX/SumdEdX;
739  dst_dedx_st dedx;
740  dedx.id_track = Id;
741  dedx.det_id = kTpcId; // TPC track
742  dedx.method = kEnsembleTruncatedMeanId; // == kTruncatedMeanId+1;
743  dedx.ndedx = TMath::Min(99,N70) + 100*((int) TrackLength);
744  dedx.dedx[0] = I70;
745  dedx.dedx[1] = D70;
746  dedx.dedx[2] = dXavLog2;
747  if ((TESTBIT(m_Mode, kCalibration))) // uncorrected dEdx
748  AddEdxTraits(tracks, dedx);
749  if (! TESTBIT(m_Mode, kDoNotCorrectdEdx)) {
750  dedx.det_id = kTpcId; // TPC track
751  m_TpcdEdxCorrection->dEdxTrackCorrection(0,dedx, etaG);
752  dedx.det_id = kTpcId; // TPC track
753  dedx.method = kTruncatedMeanId;
754  AddEdxTraits(tracks, dedx);
755  }
756  // likelihood fit
757  Double_t chisq, fitZ, fitdZ;
758 #ifdef __BENCHMARKS__DOFIT_ZN__
759  myBenchmark.Start("StdEdxY2Maker::DoFitZ");
760 #endif /* __BENCHMARKS__DOFIT_ZN__ */
761  DoFitZ(chisq, fitZ, fitdZ);
762  if (chisq >0 && chisq < 10000.0) {
763  dXavLog2 = 1;
764  SumdEdX = 0;
765  SumdX = 0;
766  for (k = 0; k < NdEdx; k++) {
767  SumdEdX += dEdxS[k].F.dEdx;
768  SumdX += dEdxS[k].F.dEdx*TMath::Log2(dEdxS[k].F.dx);
769  }
770  if (SumdEdX > 0) dXavLog2 = SumdX/SumdEdX;
771  dedx.id_track = Id;
772  dedx.det_id = kTpcId; // TPC track
773  dedx.method = kWeightedTruncatedMeanId;// == kLikelihoodFitId+1;
774  dedx.ndedx = TMath::Min(99,NdEdx) + 100*((int) TrackLength);
775  dedx.dedx[0] = TMath::Exp(fitZ);
776  dedx.dedx[1] = fitdZ;
777  dedx.dedx[2] = dXavLog2;
778  if ((TESTBIT(m_Mode, kCalibration))) // uncorrected dEdx
779  AddEdxTraits(tracks, dedx);
780  if (! TESTBIT(m_Mode, kDoNotCorrectdEdx)) {
781  dedx.det_id = kTpcId; // TPC track
782  m_TpcdEdxCorrection->dEdxTrackCorrection(2,dedx, etaG);
783  dedx.det_id = kTpcId; // TPC track
784  dedx.method = kLikelihoodFitId;
785  AddEdxTraits(tracks, dedx);
786  }
787 #ifdef __BENCHMARKS__DOFIT_ZN__
788  myBenchmark.Stop("StdEdxY2Maker::DoFitZ");
789 #endif /* __BENCHMARKS__DOFIT_ZN__ */
790 #ifdef __DEBUG1__
791 #if defined(__DEBUG_dEdx__) || defined(__DEBUG_dNdx__)
792  StThreeVectorD g3 = gTrack->geometry()->momentum(); // p of global track
793  Double_t pMomentum = g3.mag();
794  Double_t bgPion = pMomentum/StProbPidTraits::mPidParticleDefinitions[kPidPion]->mass();
795  Double_t dEdxFitL10 = TMath::Log10( 1e6*dedx.dedx[0]);
796  Double_t dNdx = StdEdxPull::EvalPred(bgPion, 2);
797  if ((TMath::Abs(pMomentum - 1.0) < 0.1 && dEdxFitL10 > 0.5 && dEdxFitL10 < 0.8) ||
798  (TMath::Abs(pMomentum - 0.5) < 0.1 && dEdxFitL10 > 1.0 && dEdxFitL10 < 1.4)) {
799  static Int_t ibreak = 0;
800  gTrack->Print();
801  if (gTrack->idTruth() > 0) {
802  St_g2t_track *g2t_track = (St_g2t_track *) GetDataSet("geant/g2t_track");
803  if (g2t_track) g2t_track->Print(gTrack->idTruth()-1,1);
804  }
805  PrintdEdx(2);
806  ibreak++;
807  }
808 #endif
809 #endif /* __DEBUG1__ */
810  if (fUsedNdx) {
811  // likelihood fit of no. of primary cluster per cm
812  Double_t chisqN, fitN = fitZ, fitdN;
813 #ifdef __BENCHMARKS__DOFIT_ZN__
814  myBenchmark.Start("StdEdxY2Maker::DoFitN");
815 #endif /* __BENCHMARKS__DOFIT_ZN__ */
816  DoFitN(chisqN, fitN, fitdN);
817  if (chisqN > -900.0 &&chisqN < 10000.0) {
818  dedx.id_track = Id;
819  dedx.det_id = kTpcId; // TPC track
820  dedx.method = kOtherMethodId2;
821  dedx.ndedx = TMath::Min(99,NdEdx) + 100*((int) TrackLength);
822  dedx.dedx[0] = fitN;
823  dedx.dedx[1] = fitdN/fitN;
824  dedx.dedx[2] = dXavLog2;
825  AddEdxTraits(tracks, dedx);
826  dedx.det_id = kTpcId; // TPC track
827  m_TpcdEdxCorrection->dEdxTrackCorrection(1,dedx, etaG);
828  dedx.det_id = kTpcId; // TPC track
829  dedx.method = kOtherMethodId;
830  AddEdxTraits(tracks, dedx);
831 #ifdef __DEBUG_dNdx__1
832  Double_t dEdxL10 = TMath::LogE()*fitZ + 6;
833  Double_t dNdxL10 = TMath::Log10(fitN);
834  if (dNdxL10 > 3.70 && dEdxL10 > -0.4) {
835  static Int_t ibreak = 0;
836  cout << "DEBUG dN/dx: dE/dx = " << TMath::Power(10.,dEdxL10) << " keV, N = " << fitN << endl;
837  PrintdEdx(2);
838  ibreak++;
839  }
840 #endif
841  }
842 #ifdef __BENCHMARKS__DOFIT_ZN__
843  myBenchmark.Stop("StdEdxY2Maker::DoFitN");
844 #endif /* __BENCHMARKS__DOFIT_ZN__ */
845  }
846  }
847 
848 #ifdef __ADD_PROB__
849  if (! TESTBIT(m_Mode, kDoNotCorrectdEdx)) {
850  StThreeVectorD g3 = gTrack->geometry()->momentum(); // p of global track
851  Double_t pMomentum = g3.mag();
852  Float_t Chisq[KPidParticles];
853  for (Int_t hyp = 0; hyp < KPidParticles; hyp++) {
854  Double_t bgL10 = TMath::Log10(pMomentum*TMath::Abs(StProbPidTraits::mPidParticleDefinitions[hyp]->charge())/StProbPidTraits::mPidParticleDefinitions[hyp]->mass());
855  Chisq[hyp] = LikeliHood(bgL10,NdEdx,FdEdx, StProbPidTraits::mPidParticleDefinitions[hyp]->charge()*StProbPidTraits::mPidParticleDefinitions[hyp]->charge());
856  }
857  for (Int_t l = 0; l < 2; l++) {if (tracks[l]) tracks[l]->addPidTraits(new StProbPidTraits(NdEdx,kTpcId,KPidParticles,Chisq));}
858  }
859 #endif /* __ADD_PROB__ */
860  }
861  } // (hvec.size() && ! TESTBIT(m_Mode, kDoNotCorrectdEdx))
862  if (pTrack) QAPlots(gTrack);
863  if ((TESTBIT(m_Mode, kCalibration))) {
864  if (! pTrack) continue; // reject non primary tracks
865 #ifdef __BEST_VERTEX__
866  if (! pEvent->primaryVertex()) continue;
867  // AuAu could have wrong ranking
868  if (pEvent->primaryVertex()->ranking() > 0) {
869  if (pTrack->vertex() != pEvent->primaryVertex()) continue; // only the first primary vertex
870  // if ( ((StPrimaryVertex *) pTrack->vertex() )->numMatchesWithBEMC() <= 0) continue;
871  } else {// try to use VpdZ to select best vertex
872  if (pVbest && pTrack->vertex() != pVbest) continue;
873  }
874 #endif /* __BEST_VERTEX__ */
875  Histogramming(gTrack);
876  }
877  }
878  if (Debug() > 1) {
879  LOG_QA << "StdEdxY2Maker:"
880  << " Type: " << pEvent->type()
881  << " Run: " << pEvent->runId()
882  << " Event: " << pEvent->id()
883  << " # track nodes: "
884  << pEvent->trackNodes().size() << endm;
885  }
886  if (mHitsUsage) mHitsUsage->Fill(TMath::Log10(TotalNoOfTpcHits+1.), TMath::Log10(NoOfTpcHitsUsed+1.));
887 #if defined(__SpaceCharge__) || defined(__CHECK_RDOMAP_AND_VOLTAGE__)
888  if ((TESTBIT(m_Mode, kCalibration))) {
889 #ifdef __SpaceCharge__
890  if (! AdcSC) {
891  AdcSC = new TH2F("AdcSC","ADC total versus z and row (-ve for East)",210,-210,210, 145, -72.5, 72.5);
892  AdcOnTrack = new TH2F("AdcOnTrack","ADC on Track versus z and row (-ve for East)",210,-210,210, 145, -72.5, 72.5);
893  dEOnTrack = new TH2F("dEOnTrack","dE (keV) on Track versus z and row (-ve for East)",210,-210,210, 145, -72.5, 72.5);
894  }
895 #endif
896 #ifdef __CHECK_RDOMAP_AND_VOLTAGE__
897  if (! AlivePads || ! ActivePads) {
898  Int_t NoOfPads = 182;
899  if (St_tpcPadConfigC::instance()->iTPC(1)) { // iTpc for all TPC sectors
900  NoOfPads = St_tpcPadConfigC::instance()->numberOfPadsAtRow(1,72);
901  }
902  Int_t nrows = St_tpcPadConfigC::instance()->numberOfRows(20);
903  if (GetTFile()) GetTFile()->cd();
904  AlivePads = new TH3F("AlivePads","Active pads from RDO map, tpcGainPadT0, and Tpc Anode Voltage:sector:row:pad",24,0.5,24.5,nrows,0.5,nrows+.5,NoOfPads,0.5,NoOfPads+0.5);
905  for (Int_t sector = 1; sector <= 24; sector++) {
906  for(Int_t row = 1; row <= St_tpcPadConfigC::instance()->numberOfRows(sector); row++) {
907  Int_t noOfPadsAtRow = St_tpcPadConfigC::instance()->numberOfPadsAtRow(sector,row);
908  if ( ! St_tpcAnodeHVavgC::instance()->livePadrow(sector,row)) continue;
909  for(Int_t pad = 1; pad<=noOfPadsAtRow; pad++) {
910  Int_t iRdo = StDetectorDbTpcRDOMasks::instance()->rdoForPadrow(sector,row,pad);
911  if ( ! StDetectorDbTpcRDOMasks::instance()->isOn(sector,iRdo)) continue;
912  Double_t gain = St_tpcPadGainT0BC::instance()->Gain(sector,row,pad);
913  if (gain <= 0.0) continue;
914  AlivePads->Fill(sector, row, pad, gain);
915  }
916  }
917  }
918  ActivePads = new TProfile3D("ActivePads","Cluster Adc:sector:row:pad",24,0.5,24.5,nrows,0.5,nrows+.5,NoOfPads,0.5,NoOfPads+0.5);
919  }
920 #endif /* __CHECK_RDOMAP_AND_VOLTAGE__ */
921  for (UInt_t i = 0; i <= TpcHitCollection->numberOfSectors(); i++) {
922  StTpcSectorHitCollection* sectorCollection = TpcHitCollection->sector(i);
923  if (sectorCollection) {
924  Int_t numberOfPadrows = sectorCollection->numberOfPadrows();
925  for (Int_t j = 0; j < numberOfPadrows; j++) {
926  StTpcPadrowHitCollection *rowCollection = sectorCollection->padrow(j);
927  if (rowCollection) {
928  StSPtrVecTpcHit &hits = rowCollection->hits();
929  Long_t NoHits = hits.size();
930  for (Long64_t k = 0; k < NoHits; k++) {
931  const StTpcHit *tpcHit = static_cast<const StTpcHit *> (hits[k]);
932  if (!tpcHit) continue;
933  Int_t sector = tpcHit->sector();
934  Int_t row = tpcHit->padrow();
935  Int_t rowc = row;
936  if (sector > 12) rowc = - row;
937  Int_t adc = tpcHit->adc();
938  Double_t Z = tpcHit->position().z();
939  AdcSC->Fill(Z,rowc, adc);
940  if (tpcHit->usedInFit()) {
941  if (tpcHit->charge() > 0) {
942  AdcOnTrack->Fill(Z,rowc, adc);
943  dEOnTrack->Fill(Z,rowc, 1e6*tpcHit->charge());
944  }
945  }
946 #ifdef __CHECK_RDOMAP_AND_VOLTAGE__
947  Int_t pad = tpcHit->pad();
948 #if 0
949  Int_t iRdo = StDetectorDbTpcRDOMasks::instance()->rdoForPadrow(sector,row,pad);
950  if ( ! StDetectorDbTpcRDOMasks::instance()->isOn(sector,iRdo)) continue;
951 #endif
952  ActivePads->Fill(sector, row, pad, adc);
953 #endif /* __CHECK_RDOMAP_AND_VOLTAGE__ */
954  }
955  }
956  }
957  }
958  }
959  }
960 #endif /* __SpaceCharge__ || __CHECK_RDOMAP_AND_VOLTAGE__ */
961 #ifdef __BENCHMARKS__DOFIT_ZN__
962  Float_t rt, cp;
963  myBenchmark.Summary(rt,cp);
964 #endif /* __BENCHMARKS__DOFIT_ZN__ */
965  return kStOK;
966 }
967 //________________________________________________________________________________
968 void StdEdxY2Maker::SortdEdx() {
969  Int_t i;
970  TArrayI idxT(NdEdx); Int_t *idx = idxT.GetArray();
971  TArrayD dT(NdEdx); Double_t *d = dT.GetArray();
972  for (i = 0; i < NdEdx; i++) d[i] = CdEdx[i].F.dEdx;
973  TMath::Sort(NdEdx,d,idx,0);
974  for (i=0;i<NdEdx;i++) dEdxS[i] = CdEdx[idx[i]];
975  TArrayI rowT(NdEdx); Int_t *r = rowT.GetArray();
976  for (i = 0; i < NdEdx; i++) r[i] = CdEdx[i].row;
977  TMath::Sort(NdEdx,r,idx,0);
978  for (i=0;i<NdEdx;i++) FdEdx[i] = CdEdx[idx[i]];
979 }
980 //________________________________________________________________________________
981 void StdEdxY2Maker::Histogramming(StGlobalTrack* gTrack) {
982  // Histograms
983  static THnSparseF *Time = 0, *TimeC = 0; // , *TimeP = 0
984  Int_t NoRows = St_tpcPadConfigC::instance()->numberOfRows(20);
985  // static Hists3D PressureT("PressureT","log(dE/dx)","row","Log(Pressure*298.2/inputGasTemperature)",-NoRows,150, 6.84, 6.99);
986 
987  // static Hists3D Volt("Volt","log(dE/dx)","Sector*Channels","Voltage", numberOfSectors*NumberOfChannels,410,990.,1400.);
988 
989  static Hists3D AvCurrent("AvCurrent","log(dEdx/Pion)","Sector*Channels","Average Current [#{mu}A]",numberOfSectors*NumberOfChannels,200,0.,1.0);
990  static Hists3D Qcm("Qcm","log(dEdx/Pion)","Sector*Channels","Accumulated Charge [uC/cm]",numberOfSectors*NumberOfChannels,200,0.,1000);
991  static Hists3D ADC3("ADC3","<logADC)>","sector","row",numberOfSectors,
992  NoRows,0,-1,
993  100,0.,10.,
994  0,-1,
995  1); // Hists3D::NtotHist = 1;
996 #ifndef MakeString
997 #define MakeString(PATH) # PATH
998 #endif
999 #ifdef __dZdY_dXdY__ /* skip dZdY and dXdY */
1000 #define __BOOK__VARS__dZdY(SIGN,NEGPOS) \
1001  static Hists3D dZdY3 ## SIGN ("dZdY3" MakeString(SIGN) ,"log(dEdx/Pion)" MakeString(NEGPOS) ,"row","dZdY",-NoRows,200,-5,5); \
1002  static Hists3D dXdY3 ## SIGN ("dXdY3" MakeString(SIGN) ,"log(dEdx/Pion)" MakeString(NEGPOS) ,"row","dXdY",-NoRows,200,-2.5,2.5);
1003 #else
1004 #define __BOOK__VARS__dZdY(SIGN,NEGPOS)
1005 #endif
1006 #ifdef __Pad_Tmbk__ /* skip Pad and Tbk */
1007 #define __BOOK__VARS__PadTmbk(SIGN,NEGPOS) \
1008  static Hists3D nPad3 ## SIGN ("nPad3" MakeString(SIGN) ,"log(dEdx/Pion)" MakeString(NEGPOS) ,"row","npad",-NoRows,18,0.5,18.5); \
1009  static Hists3D nTbk3 ## SIGN ("nTbk3" MakeString(SIGN) ,"log(dEdx/Pion)" MakeString(NEGPOS) ,"row","ntimebuckets",-NoRows,35,2.5,37.5);
1010 #else
1011 #define __BOOK__VARS__PadTmbk(SIGN,NEGPOS)
1012 #endif
1013 
1014 #define __BOOK__VARS__(SIGN,NEGPOS) \
1015  static Hists3D SecRow3 ## SIGN ("SecRow3" MakeString(SIGN) ,"<log(dEdx/Pion)>" MakeString(NEGPOS) ,"sector","row",numberOfSectors,NoRows); \
1016  static Hists3D Pressure ## SIGN ("Pressure" MakeString(SIGN) ,"log(dE/dx)" MakeString(NEGPOS) ,"row","Log(Pressure)",-NoRows,150, 6.84, 6.99); \
1017  static Hists3D Voltage ## SIGN ("Voltage" MakeString(SIGN) ,"log(dE/dx)" MakeString(NEGPOS) ,"Sector*Channels","Voltage - Voltage_{nominal}", numberOfSectors*NumberOfChannels,22,-210,10); \
1018  static Hists3D Z3 ## SIGN ("Z3" MakeString(SIGN) ,"<log(dEdx/Pion)>" MakeString(NEGPOS) ,"row","Drift Distance",-NoRows,220,-5,215); \
1019  static Hists3D G3 ## SIGN ("G3" MakeString(SIGN) ,"<log(dEdx/Pion)>" MakeString(NEGPOS) ,"row","drift time to Gating Grid (us)",-NoRows,100,-5,15); \
1020  static Hists3D xyPad3 ## SIGN ("xyPad3" MakeString(SIGN) ,"log(dEdx/Pion)" MakeString(NEGPOS) ,"sector+yrow[-0.5,0.5] and xpad [-1,1]"," xpad",numberOfSectors*20, 32,-1,1, 200, -5., 5., 0.5, 24.5); \
1021  static Hists3D dX3 ## SIGN ("dX3" MakeString(SIGN) ,"log(dEdx/Pion)" MakeString(NEGPOS) ,"row","log2(dX)",-NoRows,40,-0.5,7.5); \
1022  static Hists3D Eta3 ## SIGN ("Eta3" MakeString(SIGN) ,"log(dEdx/Pion) MC" MakeString(NEGPOS) ,"row","#eta_{G}",-NoRows,50,-2.5,2.5); \
1023  static Hists3D EtaB3 ## SIGN ("EtaB3" MakeString(SIGN) ,"log(dEdx/Pion) RC" MakeString(NEGPOS) ,"row","#eta_{G}",-NoRows,50,-2.5,2.5); \
1024 __BOOK__VARS__dZdY(SIGN,NEGPOS) \
1025 __BOOK__VARS__PadTmbk(SIGN,NEGPOS)
1026 #if 0 /* skip Pad and Tbk */
1027  static Hists3D nPad3 ## SIGN ("nPad3" MakeString(SIGN) ,"log(dEdx/Pion)" MakeString(NEGPOS) ,"row","npad",-NoRows,18,0.5,18.5); \
1028  static Hists3D nTbk3 ## SIGN ("nTbk3" MakeString(SIGN) ,"log(dEdx/Pion)" MakeString(NEGPOS) ,"row","ntimebuckets",-NoRows,35,2.5,37.5);
1029 #endif
1030  static Hists3D xyPad3qB("xyPad3qB","log(dEdx/Pion) for all","24*qB+sector+yrow[-0.5,0.5] and xpad [-1,1]"," xpad",2*numberOfSectors*20, 32,-1,1, 200, -5., 5., 0.5, 48.5);
1031 #if ! defined(__NEGATIVE_ONLY__) && ! defined(__NEGATIVE_AND_POSITIVE__)
1032  __BOOK__VARS__(,);
1033 #else
1034  __BOOK__VARS__(, for negative);
1035 #ifdef __NEGATIVE_AND_POSITIVE__
1036  __BOOK__VARS__(P, for positive);
1037 #endif
1038 #endif
1039  static TH2F *ZdcCP = 0, *BBCP = 0;
1040  // static TH2F *ctbWest = 0, *ctbEast = 0, *ctbTOFp = 0, *zdcWest = 0, *zdcEast = 0;
1041  enum {kTotalMethods = 6};
1042  // ch
1043  static TH3F *TPoints[2][kTotalMethods] = {0}; // *N[6] = {"B","70B","BU","70BU","N", "NU"}; 0 => "+", 1 => "-";
1044  static TH3F *NPoints[2][kTotalMethods] = {0}; // *N[6] = {"B","70B","BU","70BU","N", "NU"}; 0 => "+", 1 => "-";
1045  static StDedxMethod kTPoints[kTotalMethods] = {// {"F","70","FU","70U","N", "NU"};
1046  kLikelihoodFitId, // F
1047  kTruncatedMeanId, // 70
1048  kWeightedTruncatedMeanId, // FU
1049  kEnsembleTruncatedMeanId // 70U
1050  ,kOtherMethodId // N
1051  ,kOtherMethodId2 // NU
1052  };
1053 #ifdef __FIT_PULLS__
1054  static TH2F *Pulls[2][kTotalMethods] = {0};
1055  static TH2F *nPulls[2][kTotalMethods] = {0};
1056  enum {kNPulls = 3};
1057  struct PullH_t {
1058  StPidStatus::PiDStatusIDs kPid;
1059  Hists2D* histograms;
1060  };
1061  static PullH_t PullH[kNPulls] = {
1062  { StPidStatus::kI70, new Hists2D("I70")},
1063  { StPidStatus::kFit, new Hists2D("fitZ")},
1064  { StPidStatus::kdNdx, 0}
1065  };
1066  static Bool_t NotYetDone = kTRUE;
1067  if (NotYetDone) {
1068  NotYetDone = kFALSE;
1069  if (fUsedNdx) PullH[2].histograms = new Hists2D("fitN");
1070  }
1071 #endif /* __FIT_PULLS__ */
1072 #ifdef __TEST_DX__
1073  static TH3F *dXTest[2] = {0};
1074 #endif /* __TEST_DX__ */
1075  const static Int_t Nlog2dx = 80;
1076  const static Double_t log2dxLow = 0.0, log2dxHigh = 4.0;
1077  // ProbabilityPlot
1078  // static TH3F *Prob = 0;
1079  // end of ProbabilityPlot
1080  static Int_t hMade = 0;
1081  //#define __ETA_PLOTS__
1082 #ifdef __ETA_PLOTS__
1083  static TH2F *Eta[2] = {0}; // 0 -> F, 1 -> 70
1084 #endif /* __ETA_PLOTS__ */
1085  if (! gTrack && !hMade) {
1086  TFile *f = GetTFile();
1087  assert(f);
1088  f->cd();
1089  hMade=2004;
1090  // book histograms
1091  Bool_t fSetDefaultSumw2 = TH1::GetDefaultSumw2();
1092  TH1::SetDefaultSumw2(kFALSE);
1093 #ifdef __BEST_VERTEX__
1094  PVxyz = new TH3F("PVxyz","xyz for all primary vertices",100,-10,10,100,-10,10,210,-210,210);
1095  PVxyzC = new TH3F("PVxyzC","xyz for the best primary vertex",100,-10,10,100,-10,10,210,-210,210);
1096  EtaVspT[0][0] = new TH2F("EtaVspTGlP","Eta vs Log_{10}p_{T} for All positive tracks", 350, -2, 1.5, 600, -3.0, 3.0);
1097  EtaVspT[0][1] = new TH2F("EtaVspTGlN","Eta vs Log_{10}p_{T} for All negative tracks", 350, -2, 1.5, 600, -3.0, 3.0);
1098  EtaVspT[1][0] = new TH2F("EtaVspTPrP","Eta vs Log_{10}p_{T} for All positive primary tracks", 350, -2, 1.5, 600, -3.0, 3.0);
1099  EtaVspT[1][1] = new TH2F("EtaVspTPrN","Eta vs Log_{10}p_{T} for All negative primary tracks", 350, -2, 1.5, 600, -3.0, 3.0);
1100  EtaVspTC[0] = new TH2F("EtaVspTPC","Eta vs Log_{10}p_{T} for All positive global tracks used for calibration", 350, -2, 1.5, 600, -3.0, 3.0);
1101  EtaVspTC[1] = new TH2F("EtaVspTNC","Eta vs Log_{10}p_{T} for All negative global tracks used for calibration", 350, -2, 1.5, 600, -3.0, 3.0);
1102 #endif /* __BEST_VERTEX__ */
1103  mHitsUsage = new TH2F("HitsUsage","log10(No.of Used in dE/dx hits) versus log10(Total no. of Tpc Hits",
1104  80,0,8,60,0,6);
1105  Int_t nZBins = 200;
1106  Double_t ZdEdxMin = -5;
1107  Double_t ZdEdxMax = 5;
1108  ZdcCP = new TH2F("ZdcCP","ZdcCoincidenceRate (log10)",100,0,10,nZBins,ZdEdxMin,ZdEdxMax);
1109  BBCP = new TH2F("BBCP","BbcCoincidenceRate (log10)",60,0,6,nZBins,ZdEdxMin,ZdEdxMax);
1110  // TPoints block
1111  const Char_t *NS[2] = {"P",""};
1112  const Char_t *TS[2] = {"Positive","Negative"};
1113  const Char_t *N[kTotalMethods] = {"F","70","FU","70U","N", "NU"};
1114  const Char_t *T[kTotalMethods] = {"dEdx(fit)/Pion",
1115  "dEdx(I70)/Pion",
1116  "dEdx(fit_uncorrected)/Pion ",
1117  "dEdx(I70_uncorrected)/Pion",
1118  "dNdx/Pion",
1119  "dNdx(uncorrected)/Pion"};
1120  for (Int_t t = 0; t < kTotalMethods; t++) {
1121  if (! fUsedNdx && t >= 4) continue;
1122  for (Int_t s = 0; s < 2; s++) {// charge 0 => "+", 1 => "-"
1123  TPoints[s][t] = new TH3F(Form("TPoints%s%s",N[t],NS[s]),
1124  Form("%s versus Length in Tpc and <log_{2}(dX)> in TPC - iTPC %s p > 0.4 GeV/c",T[t],TS[s]),
1125  190,10,200., Nlog2dx, log2dxLow, log2dxHigh, 500,-1.,4.);
1126  NPoints[s][t] = new TH3F(Form("NPoints%s%s",N[t],NS[s]),
1127  Form("%s versus no. dEdx points in Tpc and #eta_{G} for %s p > 0.4 GeV/c",T[t],TS[s]),
1128  100,0.5,100.5, 40, -2, 2, 500,-1.,4.);
1129 #ifdef __FIT_PULLS__
1130  Pulls[s][t] = new TH2F(Form("Pull%s%s",N[t],NS[s]),
1131  Form("Pull %s versus Length in TPC %s",T[t],TS[s]),
1132  190,10.,200,nZBins,ZdEdxMin,ZdEdxMax);
1133  nPulls[s][t] = new TH2F(Form("nPull%s%s",N[t],NS[s]),
1134  Form("Pull %s versus no. dEdx points in TPC %s",T[t],TS[s]),
1135  100,0.5,100.5,nZBins,ZdEdxMin,ZdEdxMax);
1136 #endif /* __FIT_PULLS__ */
1137 #ifdef __ETA_PLOTS__
1138  if (s == 0 && t < 2) {
1139  Eta[t] = new TH2F(Form("Eta%s",N[t]),
1140  Form("%s for primary tracks versus Eta for |zPV| < 10cm and TpcLength > 40cm, TPC - iTPC",T[t]),
1141  100,-2.5,2.5,500,-1.,4.);
1142  }
1143 #endif /* __ETA_PLOTS__ */
1144  }
1145  }
1146  TDatime t1(tMin,0); // min Time and
1147  TDatime t2(tMax,0); // max
1148 
1149  UInt_t i1 = t1.Convert() - timeOffSet;
1150  UInt_t i2 = t2.Convert() - timeOffSet;
1151  Int_t Nt = (i2 - i1)/(3600); // each hour
1152  // GainMonitor = new TH2F("GainMonitor","log(dE/dx)_{corrected} - log(I(pi)) versus GainMonitor",
1153  // 100,70.,120.,nZBins,ZdEdxMin,ZdEdxMax);
1154  Int_t nBins[2] = {Nt, nZBins};
1155  Double_t xMin[2] = {(Double_t ) i1, ZdEdxMin};
1156  Double_t xMax[2] = {(Double_t ) i2, ZdEdxMax};
1157  Time = new THnSparseF("Time","log(dE/dx)_{uncorrected} - log(I(pi)) versus Date& Time", 2, nBins, xMin, xMax); f->Add(Time);
1158  TimeC = new THnSparseF("TimeC","log(dE/dx)_{corrected} - log(I(pi)) versus Date& Time after correction", 2, nBins, xMin, xMax); f->Add(TimeC);
1159  // TimeP = new THnSparseF("TimeP","log(dE/dx)_{after pressure correction} - log(I(pi)) versus Date& Time", 2, nBins, xMin, xMax); f->Add(TimeP);
1160  TH1::SetDefaultSumw2(fSetDefaultSumw2);
1161 #ifdef __TEST_DX__
1162  if (! dXTest[0]) {
1163  dXTest[0] = new TH3F("dxTestP","dX = dX_TrackFit - dX_Helix > 1e-4 versus pad row and dX_TrackFit for Positive",145,-72.5,72.5,100,-1.,9.,100,-0.25,0.25);
1164  dXTest[1] = new TH3F("dxTest" ,"dX = dX_TrackFit - dX_Helix > 1e-4 versus pad row and dX_TrackFit for Negative",145,-72.5,72.5,100,-1.,9.,100,-0.25,0.25);
1165  }
1166 #endif /* __TEST_DX__ */
1167  return;
1168  }
1169  // fill histograms
1170  tpcGas_st *tpcGas = 0;
1171  if ( m_TpcdEdxCorrection && m_TpcdEdxCorrection->tpcGas()) tpcGas = m_TpcdEdxCorrection->tpcGas()->GetTable();
1172 
1173  StThreeVectorD g3 = gTrack->geometry()->momentum(); // p of global track
1174  Double_t pMomentum = g3.mag();
1175  Double_t etaG = g3.pseudoRapidity();
1176  StPidStatus PiDs[2] = {StPidStatus(gTrack, kTRUE), StPidStatus(gTrack, kFALSE)} ;
1177  StPidStatus &PiD = fUsedx2 ? *&PiDs[0] : *&PiDs[1];
1178  if (PiD.PiDStatus < 0) return;
1179  // Double_t bg = TMath::Log10(pMomentum/StProbPidTraits::mPidParticleDefinitions[kPidPion]->mass());
1180  Int_t sCharge = 0; // positive
1181  if (gTrack->geometry()->charge() < 0) sCharge = 1; // negative
1182 #ifdef __BEST_VERTEX__
1183  EtaVspTC[sCharge]->Fill(TMath::Log10(g3.perp()), g3.pseudoRapidity());
1184 #endif /* __BEST_VERTEX__ */
1185 #ifdef __TEST_DX__
1186  if (dXTest[0]) {
1187  for (Int_t k = 0; k < NdEdx; k++) {
1188  Int_t sector = FdEdx[k].sector;
1189  Int_t row = FdEdx[k].row;
1190  Int_t rowS = row;
1191  if (sector > 12) rowS = - rowS;
1192  if (FdEdx[k].dX_TrackFit > 0 && FdEdx[k].dX_Helix > 0) {
1193  if (TMath::Abs(FdEdx[k].dX_TrackFit - FdEdx[k].dX_Helix) > 1e-4) dXTest[sCharge]->Fill(rowS, FdEdx[k].dX_TrackFit, FdEdx[k].dX_Helix - FdEdx[k].dX_TrackFit);
1194  } else {
1195  if (FdEdx[k].dX_TrackFit > 1e-4) dXTest[sCharge]->Fill(0., FdEdx[k].dX_TrackFit, 0.1);
1196  if (FdEdx[k].dX_Helix > 1e-4) dXTest[sCharge]->Fill(0., FdEdx[k].dX_Helix, -0.1);
1197  }
1198  }
1199 
1200  }
1201 #endif /* __TEST_DX__ */
1202  StDedxMethod kMethod;
1203 #ifdef __FIT_PULLS__
1204  // Pulls
1205  Int_t k = PiD.PiDkeyU3;
1206  Int_t l;
1207  for (Int_t m = 0; m < kNPulls; m++) {// I70, Ifit, dNdx
1208  if (! PullH[m].histograms) continue;
1209  if (! PiD.fStatus[PullH[m].kPid]) continue;
1210  for (l = kPidElectron; l < KPidParticles; l++) {
1211  if (PiD.fI70 && PiD.fI70->fPiD) {
1212  PullH[m].histograms->dev[l][sCharge]->Fill(PiD.bghyp[l], PiD.fStatus[PullH[m].kPid]->dev[l]);
1213  PullH[m].histograms->dev[l][ 2]->Fill(PiD.bghyp[l], PiD.fStatus[PullH[m].kPid]->dev[l]);
1214  if (k >= 0) {
1215  PullH[m].histograms->devT[l][sCharge]->Fill(PiD.bghyp[l], PiD.fStatus[PullH[m].kPid]->dev[l]);
1216  PullH[m].histograms->devT[l][ 2]->Fill(PiD.bghyp[l], PiD.fStatus[PullH[m].kPid]->dev[l]);
1217  }
1218  }
1219  }
1220  }
1221 #endif /* __FIT_PULLS__ */
1222  for (Int_t j = 0; j < kTotalMethods; j++) {
1223  kMethod = kTPoints[j];
1224  if (pMomentum > 0.4) {
1225  if (PiDs[1].dEdxStatus(kMethod)) {
1226  TPoints[sCharge][j]->Fill(PiDs[1].dEdxStatus(kMethod)->TrackLength(),PiDs[1].dEdxStatus(kMethod)->log2dX(),PiDs[1].dEdxStatus(kMethod)->dev[kPidPion]);
1227  NPoints[sCharge][j]->Fill(PiDs[1].dEdxStatus(kMethod)->N(), etaG, PiDs[1].dEdxStatus(kMethod)->dev[kPidPion]);
1228 #if 1
1229  if (Debug() > 100) {
1230  gTrack->Print();
1231  cout << "TPoints[" << sCharge << "][" << j << "] => " << TPoints[sCharge][j]->GetName() << "\t" << TPoints[sCharge][j]->GetTitle() << endl;
1232  }
1233 #endif
1234  }
1235  }
1236  if (PiD.dEdxStatus(kMethod)) {
1237 #ifdef __FIT_PULLS__
1238  Pulls[sCharge][j]->Fill(PiD.dEdxStatus(kMethod)->TrackLength(),PiD.dEdxStatus(kMethod)->devS[kPidPion]);
1239  nPulls[sCharge][j]->Fill(PiD.dEdxStatus(kMethod)->N(),PiD.dEdxStatus(kMethod)->devS[kPidPion]);
1240 #endif /* __FIT_PULLS__ */
1241 #ifdef __ETA_PLOTS__
1242  if (j < 2 && PiD.dEdxStatus(kMethod)->TrackLength() > 40) {
1243  StTrackNode *node = gTrack->node();
1244  StPrimaryTrack *pTrack = static_cast<StPrimaryTrack*>(node->track(primary));
1245  if (pTrack) {
1246  StPrimaryVertex s*primVx = (StPrimaryVertex *) pTrack->vertex();
1247  if (primVx) {
1248  if (TMath::Abs(primVx->position().z()) < 10) {
1249  StThreeVectorD P = pTrack->geometry()->helix().momentum(bField);
1250  Double_t eta = P.pseudoRapidity();
1251  Eta[j]->Fill(eta,PiD.dEdxStatus(kMethod)->dev[kPidPion]);
1252  }
1253  }
1254  }
1255  }
1256 #endif /* __ETA_PLOTS__ */
1257  }
1258  }
1259  kMethod = kLikelihoodFitId;
1260  if (PiD.dEdxStatus(kMethod) && PiD.dEdxStatus(kMethod)->TrackLength() > 20) {
1261  // if (NoFitPoints >= 20) {
1262  Int_t k;
1263  Double_t betagamma = pMomentum/StProbPidTraits::mPidParticleDefinitions[kPidPion]->mass();
1264 #if 0
1265  Double_t bgL10 = TMath::Log10(betagamma);
1266 #endif
1267  for (k = 0; k < NdEdx; k++) {
1268  Int_t sector = FdEdx[k].sector;
1269  Int_t row = FdEdx[k].row;
1270  Int_t rowS = row;
1271  if (sector > 12) rowS = - rowS;
1272 #if 0
1273  FdEdx[k].zP = // Bichsel::Instance()->GetMostProbableZ(bgL10,1.);
1274  Bichsel::Instance()->GetMostProbableZ(bgL10,TMath::Log2(FdEdx[k].F.dx)); //remove dX
1275  FdEdx[k].sigmaP = //Bichsel::Instance()->GetRmsZ(bgL10,1.);
1276  Bichsel::Instance()->GetRmsZ(bgL10,TMath::Log2(FdEdx[k].F.dx)); //remove dX
1277  Double_t predB = 1.e-6*TMath::Exp(FdEdx[k].zP);
1278  FdEdx[k].F.dEdxN = TMath::Log(FdEdx[k].F.dEdx /predB);
1279 #else
1280 #if 0
1281  if (! PiD.fdNdx) continue; //
1282  Double_t n_P = FdEdx[k].dxC*PiD.fdNdx->Pred[kPidPion];
1283  // Double_t zdEMPV = StdEdxModel::instance()->LogdEMPVGeV(n_P);//LogdEMPV(n_P) - Bichsel::Instance()->Parameterization()->MostProbableZShift();
1284 #else
1285  Double_t n_P = FdEdx[k].dxC*StdEdxModel::instance()->dNdxEff(betagamma);
1286 #endif
1287  Double_t zdEMPV = StdEdxModel::instance()->LogdEMPVGeV(n_P);//LogdEMPV(n_P) - Bichsel::Instance()->Parameterization()->MostProbableZShift();
1288  FdEdx[k].zP = zdEMPV;
1289  FdEdx[k].sigmaP = StdEdxModel::instance()->Sigma(n_P);
1290  Double_t predB = TMath::Exp(FdEdx[k].zP);
1291  FdEdx[k].F.dEdxN = TMath::Log(FdEdx[k].F.dE/predB);
1292 #endif
1293  for (Int_t l = 0; l <= StTpcdEdxCorrection::kTpcLast; l++) {
1294  if (l == StTpcdEdxCorrection::kzCorrection ||
1295  l == StTpcdEdxCorrection::kzCorrectionC) {
1296  FdEdx[k].C[l].dEdxN = FdEdx[k].F.dEdxN - (FdEdx[k].C[ StTpcdEdxCorrection::kzCorrectionC].ddEdxL +
1297  FdEdx[k].C[ StTpcdEdxCorrection::kzCorrection ].ddEdxL);
1298  } else if (l == StTpcdEdxCorrection::kTpcSecRowB ||
1299  l == StTpcdEdxCorrection::kTpcSecRowC ||
1300  l == StTpcdEdxCorrection::kTpcRowQ) {
1301  FdEdx[k].C[l].dEdxN = FdEdx[k].F.dEdxN - (FdEdx[k].C[StTpcdEdxCorrection::kTpcSecRowB].ddEdxL +
1302  FdEdx[k].C[StTpcdEdxCorrection::kTpcSecRowC].ddEdxL +
1303  FdEdx[k].C[StTpcdEdxCorrection::kTpcRowQ].ddEdxL);
1304  } else if (l == StTpcdEdxCorrection::kTpcPadMDF ||
1305  l == StTpcdEdxCorrection::kTpcPadMDC ) {
1306  FdEdx[k].C[l].dEdxN = FdEdx[k].F.dEdxN - (FdEdx[k].C[StTpcdEdxCorrection::kTpcPadMDF].ddEdxL +
1307  FdEdx[k].C[StTpcdEdxCorrection::kTpcPadMDC].ddEdxL);
1308  } else {
1309  FdEdx[k].C[l].dEdxN = FdEdx[k].F.dEdxN - FdEdx[k].C[l].ddEdxL;
1310  }
1311  if (l) FdEdx[k].C[l].dx = FdEdx[k].C[l-1].dx;
1312  }
1313  Int_t cs = NumberOfChannels*(sector-1)+FdEdx[k].channel;
1314  // if (pMomentum > pMomin && pMomentum < pMomax &&PiD.dEdxStatus(kMethod)->TrackLength() > 40 ) continue; // { // Momentum cut
1315  if (pMomentum > pMomin && pMomentum < pMomax) { // Momentum cut
1316  if (St_trigDetSumsC::instance()) {
1317  if (FdEdx[k].Zdc > 0 && ZdcCP) ZdcCP->Fill(TMath::Log10(FdEdx[k].Zdc), FdEdx[k].F.dEdxN);
1318  if (St_trigDetSumsC::instance()->bbcX() > 0) {
1319  if (BBCP) BBCP->Fill(TMath::Log10(St_trigDetSumsC::instance()->bbcX()), FdEdx[k].F.dEdxN);
1320  }
1321  }
1322  Double_t Vars[4] = {
1323  FdEdx[k].C[StTpcdEdxCorrection::kTpcSecRowB].dEdxN,
1324  FdEdx[k].F.dEdxN,
1325  0,
1326  FdEdx[k].F.dx};
1327 #if 0
1328  Double_t dEN = 0;
1329  Double_t zdEMPV = 0;
1330  if (PiD.fdNdx) {
1331  Double_t n_P = FdEdx[k].dxC*PiD.fdNdx->Pred[kPidPion];
1332  dEN = TMath::Log(FdEdx[k].F.dE); // scale to <dE/dx>_MIP = 2.4 keV/cm
1333  zdEMPV = StdEdxModel::instance()->LogdEMPV(n_P); // ? Check dx
1334  Vars[2] = dEN - zdEMPV;
1335  };
1336 #endif
1337  // SecRow3
1338  Double_t V = FdEdx[k].Voltage;
1339  Double_t VN = (row <= St_tpcPadConfigC::instance()->innerPadRows(sector)) ? V - 1170 : V - 1390;
1340  Double_t press = 0;
1341  // ADC3
1342  if (FdEdx[k].adc > 0) {
1343  Double_t ADCL = TMath::Log(FdEdx[k].adc);
1344  ADC3.Fill(sector,row,&ADCL);
1345  }
1346 
1347  if (tpcGas) {
1348  Double_t p = tpcGas->barometricPressure;
1349  if (p > 0) {
1350  press = TMath::Log(p);
1351  }
1352  Vars[0] = FdEdx[k].C[StTpcdEdxCorrection::kTpcAccumulatedQ].dEdxN;
1353  Qcm.Fill(cs,FdEdx[k].Qcm,Vars);
1354  Vars[0] = FdEdx[k].C[StTpcdEdxCorrection::kTpcCurrentCorrection].dEdxN;
1355  AvCurrent.Fill(cs,FdEdx[k].Crow,Vars);
1356  }
1357  Double_t vars[2] = {tpcTime,FdEdx[k].C[ StTpcdEdxCorrection::ktpcTime].dEdxN};
1358  if (Time) Time->Fill(vars);
1359  // if (TimeP) {vars[1] = FdEdx[k].C[StTpcdEdxCorrection::ktpcTime].dEdxN; TimeP->Fill(vars);}
1360  if (TimeC) {vars[1] = FdEdx[k].F.dEdxN; TimeC->Fill(vars);}
1361  Double_t VarsY[8] = {0};
1362 #ifdef __dZdY_dXdY__
1363 #define __FILL___VARS__dZdY(SIGN) \
1364  Vars[0] = FdEdx[k].C[StTpcdEdxCorrection::kdZdY].dEdxN; \
1365  dZdY3 ## SIGN .Fill(rowS,FdEdx[k].dZdY,Vars); \
1366  Vars[0] = FdEdx[k].C[StTpcdEdxCorrection::kdXdY].dEdxN; \
1367  dXdY3 ## SIGN .Fill(rowS,FdEdx[k].dXdY,Vars);
1368 #else
1369 #define __FILL__VARS__dZdY(SIGN)
1370 #endif
1371 #ifdef __Pad_Tmbk__ /* skip Pad and Tbk */
1372 #define __FILL__VARS__PadTmbk(SIGN) \
1373  Vars[0] = FdEdx[k].C[StTpcdEdxCorrection::knPad].dEdxN; \
1374  nPad3 ## SIGN .Fill(rowS,FdEdx[k].Npads,&Vars[1]); \
1375  Vars[0] = FdEdx[k].C[StTpcdEdxCorrection::knTbk].dEdxN; \
1376  nTbk3 ## SIGN .Fill(rowS,FdEdx[k].Ntbks,&Vars[1]);
1377 #else
1378 #define __FILL__VARS__PadTmbk(SIGN)
1379 #endif
1380 #define __FILL__VARS__(SIGN) \
1381  Vars[0] = FdEdx[k].C[StTpcdEdxCorrection::kTpcSecRowB].dEdxN; \
1382  SecRow3 ## SIGN .Fill(sector,row,Vars); \
1383  Voltage ## SIGN .Fill(cs,VN,Vars); \
1384  Vars[0] = FdEdx[k].C[StTpcdEdxCorrection::kzCorrection].dEdxN; \
1385  Z3 ## SIGN .Fill(rowS,FdEdx[k].ZdriftDistance,Vars); \
1386  Vars[0] = FdEdx[k].C[StTpcdEdxCorrection::kGatingGrid].dEdxN; \
1387  G3 ## SIGN .Fill(rowS,FdEdx[k].driftTime,Vars); \
1388  Vars[0] = FdEdx[k].C[StTpcdEdxCorrection::kTpcPadMDF].dEdxN; \
1389  xyPad3 ## SIGN .Fill(FdEdx[k].yrow,FdEdx[k].xpad, Vars); \
1390  VarsY[0] = TMath::Log2(FdEdx[k].C[StTpcdEdxCorrection::kdXCorrection].dx); \
1391  VarsY[1] = TMath::Log2(FdEdx[k].F.dx); \
1392  Vars[0] = FdEdx[k].C[StTpcdEdxCorrection::kdXCorrection].dEdxN; \
1393  dX3 ## SIGN .FillY(rowS,VarsY,Vars); \
1394  Vars[0] = FdEdx[k].C[StTpcdEdxCorrection::kEtaCorrection].dEdxN; \
1395  Eta3 ## SIGN .Fill(rowS,FdEdx[k].etaG,Vars); \
1396  Vars[0] = FdEdx[k].C[StTpcdEdxCorrection::kEtaCorrectionB].dEdxN; \
1397  EtaB3 ## SIGN .Fill(rowS,FdEdx[k].etaG,Vars); \
1398  Vars[0] = FdEdx[k].C[StTpcdEdxCorrection::ktpcPressure].dEdxN; \
1399  Pressure ## SIGN.Fill(rowS,press,Vars); \
1400 __FILL__VARS__dZdY(SIGN) \
1401 __FILL__VARS__PadTmbk(SIGN)
1402 #if ! defined(__NEGATIVE_ONLY__) && ! defined(__NEGATIVE_AND_POSITIVE__)
1403  __FILL__VARS__();
1404 #else /* ! __NEGATIVE_AND_POSITIVE__ */
1405  #if defined(__NEGATIVE_ONLY__) || defined(__NEGATIVE_AND_POSITIVE__)
1406  if (sCharge == 1) {
1407  __FILL__VARS__();
1408  }
1409  #if defined(__NEGATIVE_AND_POSITIVE__)
1410  if (sCharge == 0) {
1411  __FILL__VARS__(P);
1412  }
1413  #endif /* __NEGATIVE_AND_POSITIVE__ */
1414  #endif /* __NEGATIVE_ONLY__ || __NEGATIVE_AND_POSITIVE__ */
1415 #endif /* __NEGATIVE_AND_POSITIVE__ */
1416  Vars[0] = FdEdx[k].C[StTpcdEdxCorrection::kTpcPadMDC].dEdxN;
1417  xyPad3qB.Fill(24*FdEdx[k].qB+FdEdx[k].yrow,FdEdx[k].xpadR, Vars);
1418  } // MIP momentum cut
1419  } // loop over dEdx points
1420  }
1421  return;
1422 }
1423 //_____________________________________________________________________________
1424 void StdEdxY2Maker::PrintdEdx(Int_t iop) {
1425  const Int_t NOpts = 20;
1426  const Char_t *Names[NOpts] = {"CdEdx","FdEdx","dEdxS","dEdxU","dEdxR",
1427  "dEdxS","dEdxS","dEdxP","dEdxt","dEdxO",
1428  "dEdxM","dEdxZ","dEdxm","dEdxT","dEdxW",
1429  "dEdxC","dEdxE","dEdxp","dEdxX","dEdxd"};
1430  if (iop < 0 || iop >= NOpts) return;
1431  dEdxY2_t *pdEdx = 0;
1432  Double_t dEdx;
1433  Double_t I = 0;
1434  Int_t N70 = NdEdx - (int) (0.3*NdEdx + 0.5);
1435  Int_t N60 = NdEdx - (int) (0.4*NdEdx + 0.5);
1436  Double_t I70 = 0, I60 = 0;
1437  Double_t avrz = 0;
1438  for (Int_t i=0; i< NdEdx; i++) {
1439  dEdx = 0;
1440  if (iop == 0) {pdEdx = &CdEdx[i]; dEdx = CdEdx[i].F.dEdx;}
1441  else if (iop == 1) {pdEdx = &FdEdx[i]; dEdx = FdEdx[i].F.dEdx;}
1442  else if (iop == 2) {pdEdx = &dEdxS[i]; dEdx = dEdxS[i].F.dEdx;}
1443  else if (iop >= 3) {pdEdx = &FdEdx[i]; dEdx = FdEdx[i].C[StTpcdEdxCorrection::kUncorrected+iop-3].dEdx;}
1444  I = (i*I + 1e6*pdEdx->F.dEdx)/(i+1);
1445  // cout << Names[iop] << " " << i << " S/R " << dEdx->sector << "/" << dEdx->row
1446  // << " dEdx(keV/cm) " << 1.e6*dEdx->dEdx << " dx " << dEdx->dx
1447  // // << " dx " << dEdx->dxH
1448  // << " x[" << dEdx->xyz[0] << "," << dEdx->xyz[1] << "," << dEdx->xyz[2] << "]"
1449  // << " d[" << dEdx->xyzD[0] << "," << dEdx->xyzD[1] << "," << dEdx->xyzD[2] << "]"
1450  // << " R[" << dEdx->resXYZ[0] << "," << dEdx->resXYZ[1] << "," << dEdx->resXYZ[2] << "]"
1451  // << " Sum " << 1.e6*I << "(keV)"
1452  // << " Prob " << dEdx->Prob << endl;
1453  dEdx *= 1e6;
1454  cout << Form("%s %2i S/R %2i/%2i dEdx(keV/cm) %8.2f dx %5.2f dxC %5.2f x[%8.2f,%8.2f,%8.2f] Qcm %7.2f AvC %7.3f",
1455  Names[iop],i,pdEdx->sector,pdEdx->row,dEdx, pdEdx->F.dx ,pdEdx->dxC, pdEdx->xyz[0], pdEdx->xyz[1],
1456  pdEdx->xyz[2],pdEdx->Qcm,pdEdx->Crow);
1457  cout << Form(" d[%8.2f,%8.2f,%8.2f] Sum %8.2f Prob %8.5f", pdEdx->xyzD[0], pdEdx->xyzD[1], pdEdx->xyzD[2],
1458  I,pdEdx->Prob) << endl;
1459  if (iop == 2) {
1460  if (i < N60) I60 += dEdx;
1461  if (i < N70) I70 += dEdx;
1462  if (i == N60 - 1) {
1463  I60 /= N60;
1464  cout << " ======================= I60 \t" << I60 << endl;
1465  }
1466  if (i == N70 - 1) {
1467  I70 /= N70;
1468  cout << " ======================= I70 \t" << I70 << endl;
1469  }
1470  }
1471  avrz += TMath::Log(dEdx);
1472  }
1473  if (NdEdx) avrz /= NdEdx;
1474  cout << "mean dEdx \t" << I << "\tExp(avrz)\t" << TMath::Exp(avrz) << endl;
1475 }
1476 //________________________________________________________________________________
1477 Double_t StdEdxY2Maker::LikeliHood(Double_t Xlog10bg, Int_t NdEdx, dEdxY2_t *dEdx, Double_t chargeSq) {
1478  //SecRowMipFitpHist298P02gh1.root correction to most probable value vs log2(dx)
1479  // static const Double_t probdx2[3] = {-3.58584e-02, 4.16084e-02,-1.45163e-02};//
1480  const static Double_t ProbCut = 1.e-4;
1481  const static Double_t GeV2keV = TMath::Log(1.e-6);
1482  Double_t f = 0;
1483  for (Int_t i=0;i<NdEdx; i++) {
1484  Double_t Ylog2dx = TMath::Log2(dEdx[i].F.dx);
1485  // Double_t Ylog2dx = 1;
1486  Double_t sigmaC = 0;
1487  Double_t zMostProb = Bichsel::Instance()->GetMostProbableZ(Xlog10bg,Ylog2dx) + TMath::Log(chargeSq);
1488  Double_t sigma = Bichsel::Instance()->GetRmsZ(Xlog10bg,Ylog2dx) + sigmaC;
1489  Double_t xi = (dEdx[i].F.dEdxL - GeV2keV - zMostProb)/sigma;
1490  Double_t Phi = Bichsel::Instance()->GetProbability(Xlog10bg,Ylog2dx,xi);
1491  dEdx[i].Prob = Phi/sigma;
1492  if (dEdx[i].Prob < ProbCut) {
1493  dEdx[i].Prob = ProbCut;
1494  }
1495  f -= 2*TMath::Log( dEdx[i].Prob );
1496  }
1497  return f;
1498 }
1499 //________________________________________________________________________________
1500 void StdEdxY2Maker::Landau(Double_t x, Double_t *val){
1501  // TF1 *LandauF = new TF1("LandauF","exp([0]-0.5*((x-[1])/[2])**2+exp([3]-0.5*((x-[4])/[5])**2+exp([6]-0.5*((x-[7])/[8])**2)))",-4,10);
1502  Double_t params[9] = {
1503  -7.74975e+00,
1504  6.53414e+00,
1505  1.21524e+00,
1506  3.31409e+00,
1507  -2.58291e+00,
1508  3.51463e+00,
1509  -3.47755e+00,
1510  3.77698e-02,
1511  6.67913e-01};
1512  Double_t dev1 = (x-params[1])/params[2];
1513  Double_t dev2 = (x-params[4])/params[5];
1514  Double_t dev3 = (x-params[7])/params[8];
1515  Double_t d = TMath::Exp(params[6]-0.5*dev3*dev3);
1516  Double_t dp = -dev3/params[8]*d;
1517  Double_t c = TMath::Exp(params[3]-0.5*dev2*dev2 + d);
1518  Double_t cp = (-dev2/params[5] + dp)*c;
1519  val[0] = params[0]-0.5*dev1*dev1+c;
1520  val[1] = - dev1/params[2]+cp;
1521 }
1522 //________________________________________________________________________________
1523 void StdEdxY2Maker::fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) {
1524  Double_t Val[2];
1525  // P10
1526  // dEdxS->Draw("sigma_z:log(x)/log(2)","","prof")
1527  static Double_t sigma_p[3] = { 5.66734e-01, -1.24725e-01, 1.96085e-02};
1528  f = 0.;
1529  gin[0] = 0.;
1530  gin[1] = 0.;
1531  for (Int_t i=0;i<NdEdx; i++) {
1532  // Double_t sigma = StTpcdEdxCorrection::SumSeries(TMath::Log(FdEdx[i].dx),3,sigma_p);
1533  Double_t X = TMath::Log(FdEdx[i].F.dx);
1534  Double_t sigma = sigma_p[2];
1535  for (Int_t n = 1; n>=0; n--) sigma = X*sigma + sigma_p[n];
1536  FdEdx[i].zdev = (FdEdx[i].F.dEdxL-par[0])/sigma;
1537  Landau(FdEdx[i].zdev,Val);
1538  FdEdx[i].Prob = TMath::Exp(Val[0]);
1539  f -= Val[0];
1540  gin[0] += Val[1]/sigma;
1541  }
1542 }
1543 //________________________________________________________________________________
1544 void StdEdxY2Maker::DoFitZ(Double_t &chisq, Double_t &fitZ, Double_t &fitdZ){
1545  Double_t avz = 0;
1546  for (Int_t i=0;i<NdEdx;i++) avz += FdEdx[i].F.dEdxL;
1547  if (NdEdx>5) {
1548  avz /= NdEdx;
1549  Double_t arglist[10];
1550  Int_t ierflg = 0;
1551  m_Minuit->SetFCN(fcn);
1552  // m_Minuit->SetPrintLevel(-1);
1553  if (Debug() < 2) {
1554  arglist[0] = -1;
1555  m_Minuit->mnexcm("set print",arglist, 1, ierflg);
1556  }
1557  arglist[0] = 0.0;
1558  m_Minuit->mnexcm("set NOW",arglist, 0, ierflg);
1559  m_Minuit->mnexcm("CLEAR",arglist, 0, ierflg);
1560  arglist[0] = 0.5;
1561  m_Minuit->mnexcm("SET ERR", arglist ,1,ierflg);
1562  m_Minuit->mnparm(0, "mean", avz, 0.01,0.,0.,ierflg); //First Guess
1563  if (Debug() < 2) arglist[0] = 1.; // 1.
1564  else arglist[0] = 0.; // Check gradient
1565  m_Minuit->mnexcm("SET GRAD",arglist,1,ierflg);
1566  arglist[0] = 500;
1567  arglist[1] = 1.;
1568  m_Minuit->mnexcm("MIGRAD", arglist ,2,ierflg);
1569  m_Minuit->mnexcm("HESSE ",arglist,0,ierflg);
1570  Double_t edm,errdef;
1571  Int_t nvpar,nparx,icstat;
1572  m_Minuit->mnstat(chisq,edm,errdef,nvpar,nparx,icstat);
1573  m_Minuit->GetParameter(0, fitZ, fitdZ);
1574  }
1575  else {
1576  fitZ = fitdZ = chisq = -999.;
1577  }
1578 }
1579 //________________________________________________________________________________
1580 void StdEdxY2Maker::TrigHistos(Int_t iok) {
1581  static TProfile *BarPressure = 0, *inputTPCGasPressure = 0;
1582  static TProfile *nitrogenPressure = 0, *gasPressureDiff = 0, *inputGasTemperature = 0;
1583  static TProfile *flowRateArgon1 = 0, *flowRateArgon2 = 0, *flowRateMethane = 0;
1584  static TProfile *percentMethaneIn = 0, *ppmOxygenIn = 0, *flowRateExhaust = 0;
1585  static TProfile *ppmWaterOut = 0, *ppmOxygenOut = 0, *flowRateRecirculation = 0;
1586  // static TProfile *Center = 0, *Height = 0, *Width = 0, *CenterPressure = 0;
1587  static TH1F *Multiplicity; // mult rate
1588  static TH1F *ZdcC = 0; // ZdcCoincidenceRate
1589  static TH1F *BBC = 0; // BbcCoincidenceRate
1590  if (! iok && !BarPressure) {
1591  TDatime t1(tMin,0); // min Time and
1592  TDatime t2(tMax,0); // max
1593  UInt_t i1 = t1.Convert() - timeOffSet;
1594  UInt_t i2 = t2.Convert() - timeOffSet;
1595  Int_t Nt = (i2 - i1)/(3600); // each hour
1596  BarPressure = new TProfile("BarPressure","barometricPressure (mbar) versus time",Nt,i1,i2);
1597  inputTPCGasPressure = new TProfile("inputTPCGasPressure","inputTPCGasPressure (mbar) versus time",Nt,i1,i2);
1598  nitrogenPressure = new TProfile("nitrogenPressure","nitrogenPressure (mbar) versus time",Nt,i1,i2);
1599  gasPressureDiff = new TProfile("gasPressureDiff","gasPressureDiff (mbar) versus time",Nt,i1,i2);
1600  inputGasTemperature = new TProfile("inputGasTemperature","inputGasTemperature (degrees K) versus time",Nt,i1,i2);
1601  flowRateArgon1 = new TProfile("flowRateArgon1","flowRateArgon1 (liters/min) versus time",Nt,i1,i2);
1602  flowRateArgon2 = new TProfile("flowRateArgon2","flowRateArgon2 (liters/min) versus time",Nt,i1,i2);
1603  flowRateMethane = new TProfile("flowRateMethane","flowRateMethane (liters/min) versus time",Nt,i1,i2);
1604  percentMethaneIn = new TProfile("percentMethaneIn","percentMethaneIn (percent) versus time",Nt,i1,i2);
1605  ppmOxygenIn = new TProfile("ppmOxygenIn","ppmOxygenIn (ppm) versus time",Nt,i1,i2);
1606  flowRateExhaust = new TProfile("flowRateExhaust","flowRateExhaust (liters/min) versus time",Nt,i1,i2);
1607  ppmWaterOut = new TProfile("ppmWaterOut","ppmWaterOut (ppm) versus time",Nt,i1,i2);
1608  ppmOxygenOut = new TProfile("ppmOxygenOut","ppmOxygenOut (ppm) versus time",Nt,i1,i2);
1609  flowRateRecirculation = new TProfile("flowRateRecirculation","flowRateRecirculation (liters/min) versus time",Nt,i1,i2);
1610  // CenterPressure = new TProfile("CenterPressureP","log(center) vs log(Pressure)",150, 6.84, 6.99);
1611  // Center = new TProfile("Center","Tpc Gain Monitor center versus Time",Nt,i1,i2);
1612  // Height = new TProfile("Height","Tpc Gain Monitor height versus Time",Nt,i1,i2);
1613  // Width = new TProfile("Width","Tpc Gain Monitor width versus Time",Nt,i1,i2);
1614  // trigDetSums histograms
1615  ZdcC = new TH1F("ZdcC","ZdcCoincidenceRate (log10)",100,0,10);
1616  Multiplicity = new TH1F("Multiplicity","Multiplicity (log10)",100,0,10);
1617  BBC = new TH1F("BBC","BbcCoincidenceRate (log10)",100,0,10);
1618  } else {
1619  tpcGas_st *tpcgas = m_TpcdEdxCorrection && m_TpcdEdxCorrection->tpcGas() ? m_TpcdEdxCorrection->tpcGas()->GetTable():0;
1620  if (tpcgas) {
1621  if (BarPressure) BarPressure->Fill(tpcTime,tpcgas->barometricPressure);
1622  if (inputTPCGasPressure) inputTPCGasPressure->Fill(tpcTime,tpcgas->inputTPCGasPressure);
1623  if (nitrogenPressure) nitrogenPressure->Fill(tpcTime,tpcgas->nitrogenPressure);
1624  if (gasPressureDiff) gasPressureDiff->Fill(tpcTime,tpcgas->gasPressureDiff);
1625  if (inputGasTemperature) inputGasTemperature->Fill(tpcTime,tpcgas->inputGasTemperature);
1626  if (flowRateArgon1) flowRateArgon1->Fill(tpcTime,tpcgas->flowRateArgon1);
1627  if (flowRateArgon2) flowRateArgon2->Fill(tpcTime,tpcgas->flowRateArgon2);
1628  if (flowRateMethane) flowRateMethane->Fill(tpcTime,tpcgas->flowRateMethane);
1629  if (percentMethaneIn)
1630  percentMethaneIn->Fill(tpcTime,tpcgas->percentMethaneIn*1000./tpcgas->barometricPressure);
1631  if (ppmOxygenIn) ppmOxygenIn->Fill(tpcTime,tpcgas->ppmOxygenIn);
1632  if (flowRateExhaust) flowRateExhaust->Fill(tpcTime,tpcgas->flowRateExhaust);
1633  if (ppmWaterOut) ppmWaterOut->Fill(tpcTime,tpcgas->ppmWaterOut);
1634  if (ppmOxygenOut) ppmOxygenOut->Fill(tpcTime,tpcgas->ppmOxygenOut);
1635  if (flowRateRecirculation) flowRateRecirculation->Fill(tpcTime,tpcgas->flowRateRecirculation);
1636  }
1637  if (St_trigDetSumsC::instance()) {
1638  if (St_trigDetSumsC::instance()->zdcX() > 0 && ZdcC) ZdcC->Fill(TMath::Log10(St_trigDetSumsC::instance()->zdcX()));
1639  if (St_trigDetSumsC::instance()->bbcX() > 0 && BBC) BBC->Fill(TMath::Log10(St_trigDetSumsC::instance()->bbcX()));
1640  if (St_trigDetSumsC::instance()->mult() > 0 && Multiplicity) Multiplicity->Fill(TMath::Log10(St_trigDetSumsC::instance()->mult()));
1641  }
1642  } // (TESTBIT(m_Mode, kGASHISTOGRAMS))
1643 }
1644 //________________________________________________________________________________
1645 void StdEdxY2Maker::XyzCheck(StGlobalCoordinate *global, Int_t iokCheck) {
1646  static TH3F *XYZ = 0, *XYZbad = 0;
1647  if (! global && !XYZ) {
1648  if (Debug()) LOG_WARN << "StdEdxY2Maker::XyzCheck XYZ check Histograms" << endm;
1649  XYZ = new TH3F("XYZ","xyz for clusters",80,-200,200,80,-200,200,84,-210,210);
1650  XYZbad = new TH3F("XYZbad","xyz for clusters with mismatched sectors",
1651  80,-200,200,80,-200,200,84,-210,210);
1652  }
1653  else
1654  if (XYZ) XYZ->Fill( global->position().x(), global->position().y(), global->position().z());
1655  if (iokCheck && XYZbad) XYZbad->Fill( global->position().x(), global->position().y(), global->position().z());
1656 }
1657 //________________________________________________________________________________
1658 void StdEdxY2Maker::QAPlots(StGlobalTrack* gTrack) {
1659  static TH2F *fTdEdx[3][5];
1660  static TH2F *fqTdEdx[3];
1661  static StTpcDedxPidAlgorithm PidAlgorithm70(kTruncatedMeanId);
1662  static StTpcDedxPidAlgorithm PidAlgorithmFitZ(kLikelihoodFitId);
1663  static StTpcDedxPidAlgorithm PidAlgorithmFitN(kOtherMethodId);
1664  static StElectron* Electron = StElectron::instance();
1665  static StPionPlus* Pion = StPionPlus::instance();
1666  static StKaonPlus* Kaon = StKaonPlus::instance();
1667  static StProton* Proton = StProton::instance();
1668  static const Double_t Log10E = TMath::Log10(TMath::Exp(1.));
1669  static Int_t first=0;
1670  if (! gTrack) {
1671  TFile *f = 0;
1672  if (TESTBIT(m_Mode, kCalibration)) {
1673  f = GetTFile();
1674  if (f) f->cd();
1675  }
1676  if (!first) {
1677  fZOfBadHits = new TH1F*[fNZOfBadHits]; memset(fZOfBadHits, 0, fNZOfBadHits*sizeof(TH1F*));
1678  fZOfBadHits[0] = new TH1F("ZOfBadHits0","Total no.of rejected clusters", 100,-210,210);
1679  // AddHist(fZOfBadHits[0]);
1680  fZOfGoodHits = new TH1F("ZOfGoodHits","Z of accepted clusters",100,-210,210);
1681  fPhiOfGoodHits = new TH1F("PhiOfGoodHits","Phi of accepted clusters",100, -TMath::Pi(), TMath::Pi());
1682  fPhiOfBadHits = new TH1F("PhiOfBadHits","Phi of rejected clusters",100, -TMath::Pi(), TMath::Pi());
1683  fTracklengthInTpcTotal = new TH1F("TracklengthInTpcTotal","Total track in TPC",100,0,200);
1684  fTracklengthInTpc = new TH1F("TracklengthInTpc","Track length in TPC used for dE/dx",100,0,200);
1685  fPadTbkAll = new TH2F("PadTbkAll","no. Pads versus no. Timebuckets for all hits",35,2.5,37.,18,0.5,18.5);
1686  fPadTbkBad = new TH2F("PadTbkBad","no. Pads versus no. Timebuckets for rejected hits",35,2.5,37.,18,0.5,18.5);
1687  const Char_t *FitName[3] = {"I70","F","N"};
1688 
1689  enum {kTotalMethods = 6};
1690  for (Int_t k = 0; k < kTotalMethods/2; k++) {
1691  const Char_t *parN[5] = {"","pi","e","K","P"};
1692  const Char_t *parT[5] = {"All","|nSigmaPion| < 1","|nSigmaElectron| < 1","|nSigmaKaon| < 1","|nSigmaProton| < 1"};
1693  Int_t ny = 500;
1694  Double_t ymin = 0, ymax = 2.5;
1695  if (k == 2) {ny = 600; ymin = 0.7; ymax = 3.7;}
1696  for (Int_t t = 0; t < 5; t++) {
1697  TString Title(Form("log10(dE/dx(%s)(keV/cm)) versus log10(p(GeV/c)) for Tpc TrackLength > 40 cm %s",FitName[k],parT[t]));
1698  if (k == 2) Title = Form("log10(dN/dx) versus log10(p(GeV/c)) for Tpc TrackLength > 40 cm %s",parT[t]);
1699  fTdEdx[k][t] = new TH2F(Form("TdEdx%s%s",FitName[k],parN[t]),Title,
1700  300,-1.5,1.5, ny, ymin, ymax);
1701  fTdEdx[k][t]->SetMarkerStyle(1);
1702  fTdEdx[k][t]->SetMarkerColor(t+1);
1703  }
1704  fqTdEdx[k] = new TH2F(Form("aTdEdx%s",FitName[k]),
1705  Form("log10(dE/dx(%s)(keV/cm)) versus q*(1.5+log10(p(GeV/c))) for Tpc TrackLength > 40 cm",FitName[k]),
1706  500, -2.5, 2.5, ny, ymin, ymax);
1707  }
1708  }
1709  if (! f && !first) {
1710  // for (Int_t i = 0; i < fNZOfBadHits; i++) AddHist(fZOfBadHits[i]);
1711  AddHist(fZOfGoodHits);
1712  AddHist(fPhiOfGoodHits);
1713  AddHist(fPhiOfBadHits);
1714  AddHist(fTracklengthInTpcTotal);
1715  AddHist(fTracklengthInTpc);
1716  for (Int_t k = 0; k < 3; k++) {
1717  for (Int_t t = 0; t < 5; t++) {
1718  AddHist(fTdEdx[k][t]);
1719  }
1720  }
1721  }
1722  first = 2004;
1723  } else {
1724  StSPtrVecTrackPidTraits &traits = gTrack->pidTraits();
1725  static StDedxPidTraits *pid = 0;
1726  static Double_t TrackLength, I70, fitZ, fitN;
1727  StThreeVectorD g3 = gTrack->geometry()->momentum(); // p of global track
1728  Double_t pMomentum = g3.mag();
1729  Double_t qCharge = gTrack->geometry()->charge();
1730  Int_t k;
1731  for (UInt_t i = 0; i < traits.size(); i++) {
1732  if (! traits[i]) continue;
1733  if ( traits[i]->IsZombie()) continue;
1734  pid = dynamic_cast<StDedxPidTraits*>(traits[i]);
1735  if (pid) {
1736  if (pid->method() == kTruncatedMeanId) {
1737  I70 = pid->mean();
1738  TrackLength = pid->length();
1739  if (TrackLength < 40) continue;
1740  k = 0;
1741  fTdEdx[k][0]->Fill(TMath::Log10(pMomentum), TMath::Log10(I70)+6.);
1742  fqTdEdx[k]->Fill(qCharge*(1.5+TMath::Log10(pMomentum)), TMath::Log10(I70)+6.);
1743  const StParticleDefinition* pd = gTrack->pidTraits(PidAlgorithm70);
1744  if (pd) {
1745  if (TMath::Abs(PidAlgorithm70.numberOfSigma(Pion)) < 1) fTdEdx[k][1]->Fill(TMath::Log10(pMomentum), TMath::Log10(I70)+6.);
1746  if (TMath::Abs(PidAlgorithm70.numberOfSigma(Electron)) < 1) fTdEdx[k][2]->Fill(TMath::Log10(pMomentum), TMath::Log10(I70)+6.);
1747  if (TMath::Abs(PidAlgorithm70.numberOfSigma(Kaon)) < 1) fTdEdx[k][3]->Fill(TMath::Log10(pMomentum), TMath::Log10(I70)+6.);
1748  if (TMath::Abs(PidAlgorithm70.numberOfSigma(Proton)) < 1) fTdEdx[k][4]->Fill(TMath::Log10(pMomentum), TMath::Log10(I70)+6.);
1749  }
1750  }
1751  if (pid->method() == kLikelihoodFitId) {
1752  fitZ = TMath::Log(pid->mean()+3e-33);
1753  TrackLength = pid->length();
1754  if (TrackLength < 40) continue;
1755  k = 1;
1756  fTdEdx[k][0]->Fill(TMath::Log10(pMomentum), Log10E*fitZ + 6.);
1757  fqTdEdx[k]->Fill(qCharge*(1.5+TMath::Log10(pMomentum)), Log10E*fitZ + 6.);
1758  const StParticleDefinition* pd = gTrack->pidTraits(PidAlgorithmFitZ);
1759  if (pd) {
1760  if (TMath::Abs(PidAlgorithmFitZ.numberOfSigma(Pion)) < 1) fTdEdx[k][1]->Fill(TMath::Log10(pMomentum), Log10E*fitZ + 6.);
1761  if (TMath::Abs(PidAlgorithmFitZ.numberOfSigma(Electron)) < 1) fTdEdx[k][2]->Fill(TMath::Log10(pMomentum), Log10E*fitZ + 6.);
1762  if (TMath::Abs(PidAlgorithmFitZ.numberOfSigma(Kaon)) < 1) fTdEdx[k][3]->Fill(TMath::Log10(pMomentum), Log10E*fitZ + 6.);
1763  if (TMath::Abs(PidAlgorithmFitZ.numberOfSigma(Proton)) < 1) fTdEdx[k][4]->Fill(TMath::Log10(pMomentum), Log10E*fitZ + 6.);
1764  }
1765  }
1766  if (pid->method() == kOtherMethodId) {
1767  fitN = pid->mean();
1768  TrackLength = pid->length();
1769  if (TrackLength < 40) continue;
1770  k = 2;
1771  fTdEdx[k][0]->Fill(TMath::Log10(pMomentum),TMath::Log10(fitN));
1772  fqTdEdx[k]->Fill(qCharge*(1.5+TMath::Log10(pMomentum)),TMath::Log10(fitN));
1773  const StParticleDefinition* pd = gTrack->pidTraits(PidAlgorithmFitN);
1774  if (pd) {
1775  if (TMath::Abs(PidAlgorithmFitN.numberOfSigma(Pion)) < 1) fTdEdx[k][1]->Fill(TMath::Log10(pMomentum),TMath::Log10(fitN));
1776  if (TMath::Abs(PidAlgorithmFitN.numberOfSigma(Electron)) < 1) fTdEdx[k][2]->Fill(TMath::Log10(pMomentum),TMath::Log10(fitN));
1777  if (TMath::Abs(PidAlgorithmFitN.numberOfSigma(Kaon)) < 1) fTdEdx[k][3]->Fill(TMath::Log10(pMomentum),TMath::Log10(fitN));
1778  if (TMath::Abs(PidAlgorithmFitN.numberOfSigma(Proton)) < 1) fTdEdx[k][4]->Fill(TMath::Log10(pMomentum),TMath::Log10(fitN));
1779  }
1780  }
1781  }
1782  }
1783  }
1784 }
1785 //________________________________________________________________________________
1786 void StdEdxY2Maker::BadHit(Int_t iFlag, const StThreeVectorF &xyz) {
1787 #if 0
1788  static Int_t ibreak = 0;
1789  if (iFlag == 6 && TMath::Abs(xyz.z()) < 40) {
1790  ibreak++;
1791  }
1792 #endif
1793  static const Char_t *BadCaseses[11] =
1794  {"Total no.of rejected clusters", // 0
1795  "it is not used in track fit", // 1
1796  "it is flagged ", // 2
1797  "pad does mathc with helix prediction", // 3
1798  "dx is out interval [0.5,25]", // 4
1799  "Sector/Row gain < 0", // 6
1800  "drift distance < min || drift distance > max", // 7
1801  "dE < 0 or dx < 0", // 8
1802  "Edge effect", // 9
1803  "Anode Voltage problem" // 10
1804  };
1805  if (fZOfBadHits[0]) {
1806  fZOfBadHits[0]->Fill(xyz.z());
1807  if (fPhiOfBadHits!= 0) fPhiOfBadHits->Fill(TMath::ATan2(xyz.y(),xyz.x()));
1808  if (! fZOfBadHits[iFlag+1]) {
1809  fZOfBadHits[0]->GetDirectory()->cd();
1810  fZOfBadHits[iFlag+1] = new TH1F(*fZOfBadHits[0]);
1811  fZOfBadHits[iFlag+1]->Reset();
1812  fZOfBadHits[iFlag+1]->SetName(Form("ZOfBadHits_%i",iFlag+1));
1813  if (iFlag < 11) {
1814  fZOfBadHits[iFlag+1]->SetTitle(BadCaseses[iFlag+1]);
1815  } else {
1816  assert(m_TpcdEdxCorrection);
1817  fZOfBadHits[iFlag+1]->SetTitle(m_TpcdEdxCorrection->CorrectionStatus(iFlag-10).Title);
1818  }
1819  // AddHist(fZOfBadHits[iFlag+1]);
1820  }
1821  fZOfBadHits[iFlag+1]->Fill(xyz.z());
1822  }
1823 }
1824 //________________________________________________________________________________
1825 Int_t StdEdxY2Maker::Propagate(const StThreeVectorD &middle,const StThreeVectorD &normal,
1826  const StPhysicalHelixD &helixI, const StPhysicalHelixD &helixO,
1827  StThreeVectorD &xyz, StThreeVectorD &dirG, Double_t s[2], Double_t w[2]) {
1828  xyz = StThreeVectorD();
1829  dirG = StThreeVectorD();
1830  s[0] = helixI.pathLength(middle, normal);
1831  s[1] = helixO.pathLength(middle, normal);
1832  w[0] = w[1] = 0;
1833  Double_t sA[2] = {0};
1834  if (s[0] > 1e6 && s[1] > 1e6) {
1835  return 1;
1836  } else if (s[0] <= 1e6 && s[1] < 1e6) {
1837  sA[0] = s[0]*s[0];
1838  sA[1] = s[1]*s[1];
1839  Double_t sN = sA[0] + sA[1];
1840  w[0] = sA[0]/sN;
1841  w[1] = sA[1]/sN;
1842  } else if (s[0] <= 1e6) {
1843  w[1] = 1.;
1844  } else {
1845  w[0] = 1.;
1846  }
1847  if (w[0] > 1.e-4) {xyz += w[0]*helixO.at(s[1]); dirG += w[0]*helixO.momentumAt(s[1],bField);}
1848  if (w[1] > 1.e-4) {xyz += w[1]*helixI.at(s[0]); dirG += w[1]*helixI.momentumAt(s[0],bField);}
1849  return 0;
1850 }
1851 //________________________________________________________________________________
1852 void StdEdxY2Maker::fcnN(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) {
1853  static Int_t _debug = 0;
1854 #ifdef __DEBUG_dNdx__
1855  static TCanvas *c1 = 0;
1856  static vector<Double_t> X;
1857  static vector<Double_t> F;
1858  static vector<Double_t> E;
1859  static vector<Double_t> P;
1860  // static vector<Double_t> I;
1861  if (_debug > 0 && iflag == 1) {
1862  X.clear();
1863  F.clear();
1864  E.clear();
1865  P.clear();
1866  // I.clear();
1867  // if (!c1) c1 = new TCanvas("fcn","fcn",500,1500);
1868  if (!c1) c1 = new TCanvas("fcn","fcn",500,1000);
1869  else c1->Clear();
1870  // c1->Divide(1,3);
1871  c1->Divide(1,2);
1872  }
1873 #endif /* __DEBUG_dNdx__ */
1874  f = 0;
1875  gin[0] = 0.;
1876  Double_t dNdx = par[0]; // Mu
1877  // Double_t sigma = par[1]; // extra sigma
1878  for (Int_t i = 0; i < NdEdx; i++) {
1879  Double_t dE = FdEdx[i].F.dE;
1880  Double_t dX = FdEdx[i].dxC;
1881  Double_t Np = dNdx*dX;
1882  Double_t derivative = 0;
1883  Double_t Prob = StdEdxModel::instance()->ProbdEGeVlog(TMath::Log(dE),Np, &derivative);
1884 #ifdef __DEBUG_dNdx__
1885  if (_debug && iflag == 3) {
1886  // Double_t ee = dE + TMath::Log(1e9) -TMath::Log(Np); // to eV/Np
1887  Double_t ee = StdEdxModel::instance()->Logne(TMath::Log(dE)) -TMath::Log(Np); // to eV/Np
1888  E.push_back(ee);
1889  P.push_back(Prob);
1890  // Double_t In = StdEdxModel::instance()->GGaus()->Integral(-1.,ee);
1891  // I.push_back(In);
1892  }
1893 #endif /* __DEBUG_dNdx__ */
1894  if (Prob <= 0.0) {
1895  f += 100;
1896  FdEdx[i].Prob = 0;
1897  continue;
1898  }
1899  f -= TMath::Log(Prob);
1900  // d(dNdx) = dNp/dx
1901  gin[0] -= derivative/dX/Prob;
1902  FdEdx[i].Prob = Prob;
1903  }
1904  if (_debug > 0) {
1905  if (_debug > 2) {
1906  cout << " dNdx = " << dNdx << "\tf = " << f << endl;
1907  PrintdEdx(1);
1908  cout << "===================" << endl;
1909  }
1910 #ifdef __DEBUG_dNdx__
1911  X.push_back(dNdx);
1912  F.push_back(f);
1913  if (iflag == 3) {
1914  c1->cd(1);
1915  Int_t N = X.size();
1916  TArrayD XA(N);
1917  TArrayD YA(N);
1918  for (Int_t i = 0; i < N; i++) {
1919  XA[i] = X[i];
1920  YA[i] = F[i];
1921  }
1922  if (fdNdxGraph[0]) delete fdNdxGraph[0];
1923  fdNdxGraph[0] = new TGraph(N, XA.GetArray(), YA.GetArray());
1924  fdNdxGraph[0]->SetTitle("fcn");
1925  fdNdxGraph[0]->Draw("axp");
1926  TArrayD EA(NdEdx);
1927  TArrayD PA(NdEdx);
1928  // TArrayD IA(NdEdx);
1929  for (Int_t i = 0; i < NdEdx; i++) {
1930  EA[i] = E[i];
1931  PA[i] = P[i];
1932  // IA[i] = I[i];
1933  }
1934  if (fdNdxGraph[1]) delete fdNdxGraph[1];
1935  fdNdxGraph[1] = new TGraph(NdEdx, EA.GetArray(), PA.GetArray());
1936  fdNdxGraph[1]->SetTitle("Prob");
1937  c1->cd(2);
1938  fdNdxGraph[1]->Draw("axp");
1939  if (fdNdxGraph[2]) delete fdNdxGraph[2];
1940 // fdNdxGraph[2] = new TGraph(N, EA.GetArray(), IA.GetArray());
1941 // fdNdxGraph[2]->SetTitle("Integral");
1942 // c1->cd(3);
1943 // fdNdxGraph[2]->Draw("axp");
1944  c1->Update();
1945  }
1946 #endif /* __DEBUG_dNdx__ */
1947  }
1948 }
1949 //________________________________________________________________________________
1950 void StdEdxY2Maker::DoFitN(Double_t &chisq, Double_t &fitZ, Double_t &fitdZ){
1951  Double_t dNdx = 1e9*TMath::Exp(fitZ)/92.24/1.1; // 0.080; // 80 eV per primary interaction
1952  // Double_t dNdx = 1e9*TMath::Exp(fitZ)/80.; // 0.080; // 80 eV per primary interaction
1953  Int_t ierflg = 0;
1954  m_Minuit->SetFCN(fcnN);
1955  Double_t arglist[10] = {0};
1956  if (Debug() < 2) {
1957  arglist[0] = -1;
1958  }
1959  m_Minuit->mnexcm("set print",arglist, 1, ierflg);
1960  arglist[0] = 0.0;
1961  m_Minuit->mnexcm("set NOW",arglist, 0, ierflg);
1962  m_Minuit->mnexcm("CLEAR",arglist, 0, ierflg);
1963  arglist[0] = 0.5;
1964  m_Minuit->mnexcm("SET ERR", arglist ,1,ierflg);
1965  // m_Minuit->mnparm(0, "LogdNdx", TMath::Log(dNdx), 0.5, 0.,0.,ierflg); //First Guess
1966  m_Minuit->DefineParameter(0, "dNdx", dNdx, 0.5, 0.2*dNdx, 5*dNdx);
1967  // m_Minuit->DefineParameter(1, "sigma", 0.01, 0.01, 0.0, 0.5);
1968  arglist[0] = 1.0;
1969  m_Minuit->mnexcm("CALLfcn", arglist ,1,ierflg);
1970  if (Debug() < 4) arglist[0] = 1.; // 1.
1971  else arglist[0] = 0.; // Check gradient
1972  m_Minuit->mnexcm("SET GRAD",arglist,1,ierflg);
1973  arglist[0] = 500;
1974  arglist[1] = 1.;
1975  // m_Minuit->mnexcm("MIGRAD", arglist ,2,ierflg);
1976  m_Minuit->mnexcm("MINIMIZE", arglist ,2,ierflg);
1977  m_Minuit->mnexcm("HESSE ",arglist,0,ierflg);
1978  arglist[0] = 3.0;
1979  m_Minuit->mnexcm("CALLfcn", arglist ,1,ierflg);
1980  Double_t edm,errdef;
1981  Int_t nvpar,nparx,icstat;
1982  m_Minuit->mnstat(chisq,edm,errdef,nvpar,nparx,icstat);
1983  m_Minuit->GetParameter(0, fitZ, fitdZ);
1984 }
1985 //________________________________________________________________________________
1986 void StdEdxY2Maker::IntegrateAdc(const StTpcHitCollection* TpcHitCollection) {
1987  if (! fIntegratedAdc) fIntegratedAdc = new TH2F("AdcI","Integrated Adc for timebuckets foreach High Anode Vltage socket",451,-0.5,450.5,192,0.5,192.5);
1988  else fIntegratedAdc->Reset();
1989 
1990  if (! TpcHitCollection) { cout << "No TPC Hit Collection" << endl; return;}
1991  UInt_t numberOfSectors = TpcHitCollection->numberOfSectors();
1992  for (UInt_t i = 0; i< numberOfSectors; i++) {
1993  const StTpcSectorHitCollection* sectorCollection = TpcHitCollection->sector(i);
1994  if (sectorCollection) {
1995  Int_t sector = i;
1996  Int_t numberOfPadrows = sectorCollection->numberOfPadrows();
1997  for (int j = 0; j< numberOfPadrows; j++) {
1998  const StTpcPadrowHitCollection *rowCollection = sectorCollection->padrow(j);
1999  if (rowCollection) {
2000  Int_t row = j + 1;
2001  Int_t channel = St_TpcAvgPowerSupplyC::instance()->ChannelFromRow(sector,row);
2002  Double_t sc = NumberOfChannels*(sector-1) + channel;
2003  const StSPtrVecTpcHit &hits = rowCollection->hits();
2004  Long64_t NoHits = hits.size();
2005  for (Long64_t k = 0; k < NoHits; k++) {
2006  const StTpcHit *tpcHit = static_cast<const StTpcHit *> (hits[k]);
2007  if (! tpcHit) continue;
2008  fIntegratedAdc->Fill(tpcHit->timeBucket(), sc, tpcHit->adc());
2009  }
2010  }
2011  }
2012  }
2013  }
2014  if (fIntegratedAdc->GetEntries()) {
2015  Int_t ny = fIntegratedAdc->GetNbinsY();
2016  Int_t nx = fIntegratedAdc->GetNbinsX();
2017  for (Int_t iy = 1; iy <= ny; iy++) {
2018  Double_t sum = 0;
2019  for (Int_t ix = 1; ix <= nx; ix++) {
2020  sum += fIntegratedAdc->GetBinContent(ix,iy);
2021  if (sum > 0) fIntegratedAdc->SetBinContent(ix,iy, sum);
2022  }
2023  }
2024  }
2025 }
2026 //________________________________________________________________________________
2027 Double_t StdEdxY2Maker::IntegratedAdc(const StTpcHit* tpcHit) {
2028  if (! fIntegratedAdc) return 0;
2029  Int_t sector = tpcHit->sector();
2030  Int_t row = tpcHit->padrow();
2031  Int_t channel = St_TpcAvgPowerSupplyC::instance()->ChannelFromRow(sector,row);
2032  Double_t sc = NumberOfChannels*(sector-1) + channel;
2033  return fIntegratedAdc->Interpolate(tpcHit->timeBucket(), sc);
2034 }
StMaker::m_Mode
Int_t m_Mode
counters
Definition: StMaker.h:81
tof
Definition: tof.h:15
StdEdxY2Maker::Make
virtual Int_t Make()
Definition: StdEdxY2Maker.cxx:280
StPidStatus::PiDkeyU3
Int_t PiDkeyU3
only one with devZs&lt;3,
Definition: StPidStatus.h:134
StHelix::pathLength
pair< double, double > pathLength(double r) const
path length at given r (cylindrical r)
Definition: StHelix.cc:351
dEdxY2_t
Definition: StTpcdEdxCorrection.h:141
dst_dedx_st
Definition: StTrackPidTraits.h:56
kStOK
Definition: Stypes.h:40
StMaker::Finish
virtual Int_t Finish()
Definition: StMaker.cxx:776
StdEdxY2Maker::Finish
virtual Int_t Finish()
Definition: StdEdxY2Maker.cxx:263
Generated by   doxygen 1.8.5