00001
00002 #define CompareWithToF
00003 #include <Stiostream.h>
00004 #include "StdEdxY2Maker.h"
00005
00006 #include "TMinuit.h"
00007 #include "TSystem.h"
00008 #include "TMath.h"
00009 #include "TH2.h"
00010 #include "TH3.h"
00011 #include "TStyle.h"
00012 #include "TProfile.h"
00013 #include "TProfile2D.h"
00014 #include "TTree.h"
00015 #include "TChain.h"
00016 #include "TFile.h"
00017 #include "TNtuple.h"
00018 #include "TCanvas.h"
00019 #include "TRandom.h"
00020 #include "TClonesArray.h"
00021 #include "TArrayI.h"
00022 #include "TArrayD.h"
00023
00024 #include "StMagF.h"
00025 #include "StDbUtilities/StMagUtilities.h"
00026 #include "StMessMgr.h"
00027 #include "BetheBloch.h"
00028 #include "StBichsel/Bichsel.h"
00029 #include "StDetectorId.h"
00030 #include "StDedxMethod.h"
00031
00032 #include "StTimer.hh"
00033 #include "SystemOfUnits.h"
00034 #ifndef ST_NO_NAMESPACES
00035 using namespace units;
00036 #endif
00037 #include "StarClassLibrary/StParticleTypes.hh"
00038
00039
00040 #include "StDbUtilities/StTpcCoordinateTransform.hh"
00041 #include "StDbUtilities/StCoordinates.hh"
00042 #include "StTpcDb/StTpcDb.h"
00043
00044 #include "StEventTypes.h"
00045 #include "StProbPidTraits.h"
00046 #ifdef CompareWithToF
00047 #include "StBTofPidTraits.h"
00048 #endif
00049 #include "StTpcDedxPidAlgorithm.h"
00050 #include "StDetectorDbMaker/St_tss_tssparC.h"
00051 #include "StDetectorDbMaker/St_tpcAnodeHVavgC.h"
00052 #include "StDetectorDbMaker/St_tpcAvCurrentC.h"
00053 #include "StDetectorDbMaker/St_TpcAvgCurrentC.h"
00054 #include "StDetectorDbMaker/St_trigDetSumsC.h"
00055 const static StPidParticle NHYPS = KPidParticles;
00056 const static Int_t tZero= 19950101;
00057 const static Int_t tMin = 20090301;
00058 const static Int_t tMax = 20120705;
00059 const static TDatime t0(tZero,0);
00060 const static Int_t timeOffSet = t0.Convert();
00061 const static Int_t NdEdxMax = 60;
00062 Int_t StdEdxY2Maker::NdEdx = 0;
00063 dEdxY2_t *StdEdxY2Maker::CdEdx = 0;
00064 dEdxY2_t *StdEdxY2Maker::FdEdx = 0;
00065 dEdxY2_t *StdEdxY2Maker::dEdxS = 0;
00066 static Int_t numberOfSectors = 0;
00067 static Int_t numberOfTimeBins = 0;
00068 static Int_t NumberOfRows = 0;
00069 static Int_t NumberOfInnerRows = 0;
00070 static Int_t NoPads = 0;
00071 static Double_t innerSectorPadPitch = 0;
00072 static Double_t outerSectorPadPitch = 0;
00073
00074 const static Double_t pMomin = 0.4;
00075 const static Double_t pMomax = 0.5;
00076 Bichsel *StdEdxY2Maker::m_Bichsel = 0;
00077 #include "dEdxTrackY2.h"
00078
00079
00080 const static Int_t fNZOfBadHits = 11;
00081 static TH1F **fZOfBadHits = 0;
00082 static TH1F *fZOfGoodHits = 0;
00083 static TH1F *fPhiOfBadHits = 0;
00084 static TH1F *fTracklengthInTpcTotal = 0;
00085 static TH1F *fTracklengthInTpc = 0;
00086
00087 ClassImp(StdEdxY2Maker);
00088
00089 StdEdxY2Maker::StdEdxY2Maker(const char *name):
00090 StMaker(name),
00091 m_Minuit(0), m_TpcdEdxCorrection(0), m_Mask(-1),
00092 mHitsUsage(0)
00093 {
00094 memset (beg, 0, end-beg);
00095
00096 SETBIT(m_Mode,kPadSelection);
00097 m_Minuit = new TMinuit(2);
00098 memset(mNormal[0] ,0,sizeof(mNormal ));
00099 memset(mRowPosition[0][0],0,sizeof(mRowPosition));
00100 memset(mPromptNormal[0] ,0,sizeof(mPromptNormal ));
00101 memset(mPromptPosition[0][0],0,sizeof(mPromptPosition));
00102 }
00103
00104 Int_t StdEdxY2Maker::Init(){
00105 if (IAttr("EmbeddingShortCut")) {
00106 m_Mode = 0;
00107 SETBIT(m_Mode,kEmbedding);
00108 SETBIT(m_Mode,kPadSelection);
00109 SETBIT(m_Mask,StTpcdEdxCorrection::kTpcLast);
00110 }
00111 if (Debug()) {
00112 if (! TESTBIT(m_Mode, kOldClusterFinder))
00113 gMessMgr->Warning() << "StdEdxY2Maker::Init use new Cluster Finder parameterization" << endm;
00114 else
00115 gMessMgr->Warning() << "StdEdxY2Maker::Init use old Cluster Finder parameterization" << endm;
00116 if (TESTBIT(m_Mode, kPadSelection))
00117 gMessMgr->Warning() << "StdEdxY2Maker::Init Pad Selection is ON" << endm;
00118 if (TESTBIT(m_Mode, kDoNotCorrectdEdx))
00119 gMessMgr->Warning() << "StdEdxY2Maker::Init Don't Correct dEdx" << endm;
00120 if (TESTBIT(m_Mode, kEmbedding))
00121 gMessMgr->Warning() << "StdEdxY2Maker::Init This is embedding run" << endm;
00122 }
00123 if (! m_Bichsel) m_Bichsel = Bichsel::Instance();
00124
00125 gMessMgr->SetLimit("StdEdxY2Maker:: mismatched Sector",20);
00126 gMessMgr->SetLimit("StdEdxY2Maker:: pad/TimeBucket out of range:",20);
00127 gMessMgr->SetLimit("StdEdxY2Maker:: Helix Prediction",20);
00128 gMessMgr->SetLimit("StdEdxY2Maker:: Coordinates",20);
00129 gMessMgr->SetLimit("StdEdxY2Maker:: Prediction",20);
00130 gMessMgr->SetLimit("StdEdxY2Maker:: NdEdx",20);
00131 gMessMgr->SetLimit("StdEdxY2Maker:: Illegal time for scalers",20);
00132 return StMaker::Init();
00133 }
00134
00135 Int_t StdEdxY2Maker::InitRun(Int_t RunNumber){
00136 static Int_t DoOnce = 0;
00137 if (!gStTpcDb) {
00138 cout << "Database Missing! Can't initialize StdEdxY2Maker" << endl;
00139 return kStFatal;
00140 }
00141
00142 numberOfSectors = gStTpcDb->Dimensions()->numberOfSectors();
00143 numberOfTimeBins = gStTpcDb->Electronics()->numberOfTimeBins();
00144 NumberOfRows = gStTpcDb->PadPlaneGeometry()->numberOfRows();
00145 NumberOfInnerRows = gStTpcDb->PadPlaneGeometry()->numberOfInnerRows();
00146 NoPads = TMath::Max(gStTpcDb->PadPlaneGeometry()->numberOfPadsAtRow(NumberOfInnerRows),
00147 gStTpcDb->PadPlaneGeometry()->numberOfPadsAtRow(NumberOfRows));
00148 innerSectorPadPitch = gStTpcDb->PadPlaneGeometry()->innerSectorPadPitch();
00149 outerSectorPadPitch = gStTpcDb->PadPlaneGeometry()->outerSectorPadPitch();
00150
00151 if ( ! St_trigDetSumsC::instance() ) gMessMgr->Error() << "StdEdxY2Maker:: Cannot find trigDetSums" << endm;
00152 else {
00153 if (!St_trigDetSumsC::instance()->GetNRows()) gMessMgr->Error() << "StdEdxY2Maker:: trigDetSums has not data" << endm;
00154 else {
00155 UInt_t date = GetDateTime().Convert();
00156 if (date < St_trigDetSumsC::instance()->timeOffset()) {
00157 gMessMgr->Error() << "StdEdxY2Maker:: Illegal time for scalers = "
00158 << St_trigDetSumsC::instance()->timeOffset() << "/" << date
00159 << " Run " << St_trigDetSumsC::instance()->runNumber() << "/" << GetRunNumber() << endm;
00160 }
00161 }
00162 }
00163 if (! DoOnce) {
00164 DoOnce = 1;
00165 if (TESTBIT(m_Mode, kCalibration)) {
00166 if (Debug()) gMessMgr->Warning() << "StdEdxY2Maker::InitRun Calibration Mode is On (make calibration histograms)" << endm;
00167 TFile *f = GetTFile();
00168 if (f) {
00169 f->cd();
00170 if ((TESTBIT(m_Mode, kGASHISTOGRAMS))) {
00171 if (Debug()) gMessMgr->Warning() << "StdEdxY2Maker::InitRun Gas Histograms is ON" << endm;
00172 TrigHistos();
00173 }
00174 Histogramming();
00175 if ((TESTBIT(m_Mode, kXYZcheck))) XyzCheck();
00176 }
00177 }
00178 QAPlots(0);
00179 }
00180 SafeDelete(m_TpcdEdxCorrection);
00181 m_TpcdEdxCorrection = new StTpcdEdxCorrection(m_Mask, Debug());
00182
00183 StTpcCoordinateTransform transform(gStTpcDb);
00184 for (Int_t sector = 1; sector<= numberOfSectors; sector++) {
00185 for (Int_t row = 1; row <= NumberOfRows; row++) {
00186
00187 if (Debug()>1) cout << "========= sector/row ========" << sector << "/" << row << endl;
00188 StTpcLocalSectorDirection dirLS(0.,1.,0.,sector,row); if (Debug()>1) cout << "dirLS\t" << dirLS << endl;
00189 StTpcLocalDirection dirL;
00190 StTpcLocalSectorAlignedDirection dirLSA;
00191 transform(dirLS,dirLSA); if (Debug()>1) cout << "dirLSA\t" << dirLSA << endl;
00192 transform(dirLSA,dirL); if (Debug()>1) cout << "dirL\t" << dirL << endl;
00193 StGlobalDirection dirG;
00194 transform(dirL,dirG); if (Debug()>1) cout << "dirG\t" << dirG << endl;
00195 SafeDelete(mNormal[sector-1][row-1]);
00196 mNormal[sector-1][row-1] = new StThreeVectorD(dirG.position().unit());
00197 if (Debug()>1) cout << "Normal[" << sector-1 << "][" << row-1 << "] = " << *mNormal[sector-1][row-1] << endl;
00198 Double_t padlength;
00199 if (row<14) padlength = gStTpcDb->PadPlaneGeometry()->innerSectorPadLength();
00200 else padlength = gStTpcDb->PadPlaneGeometry()->outerSectorPadLength();
00201 for (Int_t l = 0; l < 3; l++) {
00202 SafeDelete(mRowPosition[sector-1][row-1][l]);
00203 Double_t y = transform.yFromRow(row);
00204 if (l == 1) y += padlength/2.;
00205 if (l == 2) y -= padlength/2.;
00206 StTpcLocalSectorCoordinate lsCoord(0., y, 10.,sector,row); if (Debug()>1) cout << lsCoord << endl;
00207 StGlobalCoordinate gCoord;
00208 StTpcLocalSectorAlignedCoordinate lsCoordA;
00209 transform(lsCoord,lsCoordA); if (Debug()>1) cout << lsCoordA << endl;
00210 transform(lsCoordA, gCoord); if (Debug()>1) cout << gCoord << endl;
00211 SafeDelete(mRowPosition[sector-1][row-1][l]);
00212 mRowPosition[sector-1][row-1][l] =
00213 new StThreeVectorD(gCoord.position().x(),gCoord.position().y(),gCoord.position().z());
00214 if (Debug()>1) cout << "mRowPosition[" << sector-1 << "][" << row-1 << "][" << l << "] = "
00215 << *mRowPosition[sector-1][row-1][l] << endl;
00216 }
00217 }
00218 }
00219 for (Int_t iWestEast = 0; iWestEast < 2; iWestEast++) {
00220 for (Int_t io = 0; io < 2; io++) {
00221 if (Debug()>1) cout << "========= West (0) or East(1) / Inner(0) or Outer (1) ========"
00222 << iWestEast << "/" << io << endl;
00223 Int_t sector = (iWestEast == 0) ? 12 : 24;
00224 Int_t row = (io == 0) ? 1 : 14;
00225 StTpcLocalSectorDirection dirLS(0.,0.,1.0,sector,row); if (Debug()>1) cout << "dirLS\t" << dirLS << endl;
00226 StTpcLocalDirection dirL;
00227 StTpcLocalSectorAlignedDirection dirLSA;
00228 transform(dirLS,dirLSA); if (Debug()>1) cout << "dirLSA\t" << dirLSA << endl;
00229 transform(dirLSA,dirL); if (Debug()>1) cout << "dirL\t" << dirL << endl;
00230 StGlobalDirection dirG;
00231 transform(dirL,dirG); if (Debug()>1) cout << "dirG\t" << dirG << endl;
00232 SafeDelete(mPromptNormal[iWestEast][io]);
00233 mPromptNormal[iWestEast][io] = new StThreeVectorD(dirG.position().unit());
00234 if (Debug()>1) cout << "mPromptNormal[" << iWestEast << "][" << io << "] = " << *mPromptNormal[iWestEast][io] << endl;
00235
00236 Double_t z[2][3] = {
00237
00238 { -0.6 - 0.2, 0, -0.6 - 2*0.2},
00239 { -0.6 - 0.4, 0, -0.6 - 2*.04}
00240 };
00241 for (Int_t l = 0; l < 3; l++) {
00242 SafeDelete(mPromptPosition[iWestEast][io][l]);
00243 Double_t y = transform.yFromRow(row);
00244 StTpcLocalSectorCoordinate lsCoord(0., y, z[io][l],sector,row); if (Debug()>1) cout << lsCoord << endl;
00245 StGlobalCoordinate gCoord;
00246 StTpcLocalSectorAlignedCoordinate lsCoordA;
00247 transform(lsCoord,lsCoordA); if (Debug()>1) cout << lsCoordA << endl;
00248 transform(lsCoordA, gCoord); if (Debug()>1) cout << gCoord << endl;
00249 SafeDelete(mPromptPosition[iWestEast][io][l]);
00250 mPromptPosition[iWestEast][io][l] =
00251 new StThreeVectorD(gCoord.position().x(),gCoord.position().y(),gCoord.position().z());
00252 if (Debug()>1) cout << "mPromptPosition[" << sector-1 << "][" << row-1 << "][" << l << "] = "
00253 << *mPromptPosition[iWestEast][io][l] << endl;
00254 }
00255 }
00256 }
00257 return kStOK;
00258 }
00259
00260 Int_t StdEdxY2Maker::FinishRun(Int_t OldRunNumber) {
00261
00262
00263 for (int i = 0; i < 24; i++)
00264 for (int j = 0; j < 45; j++) {
00265 SafeDelete(mNormal[i][j]);
00266 for (Int_t k = 0; k < 3; k++)
00267 SafeDelete(mRowPosition[i][j][k]);
00268 }
00269
00270 SafeDelete(m_TpcdEdxCorrection);
00271 return StMaker::FinishRun(OldRunNumber);
00272 }
00273
00274 Int_t StdEdxY2Maker::Finish() {
00275 FinishRun(0);
00276 SafeDelete(m_TpcdEdxCorrection);
00277 SafeDelete(m_Minuit);
00278 SafeDelete(m_Bichsel);
00279 return StMaker::Finish();
00280 }
00281
00282 Int_t StdEdxY2Maker::Make(){
00283 static const Double_t RA[2] = { 154.484, 81.42};
00284 static const Double_t WireLenth[2] = { 3.6e5, 1.6e5};
00285 static StTimer timer;
00286 static StTpcLocalSectorCoordinate localSect[4];
00287 static StTpcPadCoordinate PadOfTrack, Pad;
00288 static StTpcLocalSectorDirection localDirectionOfTrack;
00289 static StTpcLocalSectorAlignedDirection localADirectionOfTrack;
00290 static StTpcLocalSectorAlignedCoordinate localA, lANext;
00291 static StThreeVectorD xyz[4];
00292 static StThreeVectorD dirG;
00293 static Double_t s[2], s_in[2], s_out[2], w[2], w_in[2], w_out[2], s_inP[2], s_outP[2], dx;
00294 if (Debug() > 0) timer.start();
00295 static dEdxY2_t CdEdxT[180];
00296 CdEdx = CdEdxT;
00297 FdEdx = CdEdxT + NdEdxMax;
00298 dEdxS = CdEdxT + 2*NdEdxMax;
00299 St_tpcGas *tpcGas = m_TpcdEdxCorrection->tpcGas();
00300 if (TESTBIT(m_Mode, kCalibration) && tpcGas) TrigHistos(1);
00301 StTpcCoordinateTransform transform(gStTpcDb);
00302 Double_t bField = 0;
00303 StEvent* pEvent = dynamic_cast<StEvent*> (GetInputDS("StEvent"));
00304 if (!pEvent) {
00305 gMessMgr->Info() << "StdEdxY2Maker: no StEvent " << endm;
00306 return kStOK;
00307 }
00308 if (pEvent->runInfo()) bField = pEvent->runInfo()->magneticField()*kilogauss;
00309
00310 if (fabs(bField) < 1.e-5*kilogauss) return kStOK;
00311
00312
00313 Int_t TotalNoOfTpcHits = 0;
00314 Int_t NoOfTpcHitsUsed = 0;
00315 StTpcHitCollection* TpcHitCollection = pEvent->tpcHitCollection();
00316 if (! TpcHitCollection) {
00317 gMessMgr->Info() << "StdEdxY2Maker: no TpcHitCollection " << endm;
00318 return kStOK;
00319 }
00320 TotalNoOfTpcHits = TpcHitCollection->numberOfHits();
00321 StSPtrVecTrackNode& trackNode = pEvent->trackNodes();
00322 UInt_t nTracks = trackNode.size();
00323 StTrackNode *node=0;
00324 for (UInt_t i=0; i < nTracks; i++) {
00325 node = trackNode[i];
00326 if (!node) continue;
00327 StGlobalTrack *gTrack = static_cast<StGlobalTrack *>(node->track(global));
00328 StPrimaryTrack *pTrack = static_cast<StPrimaryTrack*>(node->track(primary));
00329
00330 #if 1
00331 if (TESTBIT(m_Mode, kCalibration)) {
00332 if (! pTrack) continue;
00333 #if 1
00334 if (pTrack->vertex() != pEvent->primaryVertex()) continue;
00335 #else
00336 if ( ((StPrimaryVertex *) pTrack->vertex() )->numMatchesWithBEMC() <= 0) continue;
00337 #endif
00338 }
00339 #endif
00340
00341 StTrack *track = 0;
00342 StTrack *tracks[2] = {gTrack, pTrack};
00343 if (! TESTBIT(m_Mode, kDoNotCorrectdEdx)) {
00344
00345
00346 for (int l = 0; l < 2; l++) {
00347 track = tracks[l];
00348 if (track) {
00349 StSPtrVecTrackPidTraits &traits = track->pidTraits();
00350 UInt_t size = traits.size();
00351 if (size) {
00352 for (UInt_t i = 0; i < size; i++) {
00353 StDedxPidTraits *pid = dynamic_cast<StDedxPidTraits*>(traits[i]);
00354 if (! pid) continue;
00355 if (pid->detector() != kTpcId) continue;
00356 traits[i]->makeZombie(1);
00357 }
00358 }
00359 }
00360 }
00361 }
00362 if (! gTrack || gTrack->flag() <= 0) continue;
00363 StPhysicalHelixD helixO = gTrack->outerGeometry()->helix();
00364 StPhysicalHelixD helixI = gTrack->geometry()->helix();
00365 if (Debug() > 1) {
00366 cout << "Track:" << i
00367 << "\ttype " << gTrack->type()
00368 << "\tvertex " << gTrack->vertex()
00369 << "\tkey " << gTrack->key()
00370 << "\tflag " << gTrack->flag()
00371 << "\tencodedMethod " << gTrack->encodedMethod()
00372 << "\timpactParameter " << gTrack->impactParameter()
00373 << "\tlength " << gTrack->length()
00374 << "\tnumberOfPossiblePoints " << gTrack->numberOfPossiblePoints() << endl;
00375 cout << "pxyzI: " << gTrack->geometry()->momentum() << "\tmag " << gTrack->geometry()->momentum().mag() << endl;
00376 cout << "pxyzO: " << gTrack->outerGeometry()->momentum() << "\tmag " << gTrack->outerGeometry()->momentum().mag() << endl;
00377 cout << "start Point: " << helixI.at(0) << endl;
00378 cout << "end Point: " << helixO.at(0) << endl;
00379 }
00380
00381 StPtrVecHit hvec = gTrack->detectorInfo()->hits();
00382
00383 if (hvec.size() && ! TESTBIT(m_Mode, kDoNotCorrectdEdx)) {
00384 Int_t Id = gTrack->key();
00385 Int_t NoFitPoints = gTrack->fitTraits().numberOfFitPoints(kTpcId);
00386 NdEdx = 0;
00387 Double_t TrackLength70 = 0;
00388 Double_t TrackLength = 0;
00389 Double_t TrackLengthTotal = 0;
00390 for (UInt_t j=0; j<hvec.size(); j++) {
00391 if (hvec[j]->detector() != kTpcId) continue;
00392 StTpcHit *tpcHit = static_cast<StTpcHit *> (hvec[j]);
00393 if (! tpcHit) continue;
00394 if (Debug() > 1) {tpcHit->Print();}
00395 Int_t sector = tpcHit->sector();
00396 Int_t row = tpcHit->padrow();
00397 if (! St_tpcAnodeHVavgC::instance()->livePadrow(sector,row)) continue;
00398 xyz[3] = StThreeVectorD(tpcHit->position().x(),tpcHit->position().y(),tpcHit->position().z());
00399
00400 const StThreeVectorD &normal = *mNormal[sector-1][row-1];
00401 const StThreeVectorD &middle = *mRowPosition[sector-1][row-1][0];
00402 const StThreeVectorD &upper = *mRowPosition[sector-1][row-1][1];
00403 const StThreeVectorD &lower = *mRowPosition[sector-1][row-1][2];
00404
00405 Double_t V = St_tpcAnodeHVavgC::instance()->voltagePadrow(sector,row);
00406 if ((row <= 13 && 1170 - V > 100) ||
00407 (row > 13 && 1390 - V > 100)) {BadHit(9,tpcHit->position()); continue;}
00408
00409 if (Propagate(middle,normal,helixI,helixO,bField,xyz[0],dirG,s,w)) {BadHit(2,tpcHit->position()); continue;}
00410 if (Debug() > 1) {
00411 cout << " Prediction:\t" << xyz[0]
00412 << "\tat s=\t" << s[0] << "/" << s[1]
00413 << "\tw = " << w[0] << "/" << w[1] << endl;
00414 }
00415 StThreeVectorD dif = xyz[3] - xyz[0];
00416 if (dif.perp() > 2.0) {if (Debug() > 1) {cout << "Prediction is to far from hit:\t" << xyz[3] << endl;}
00417 continue;
00418 }
00419 if (Propagate(upper,normal,helixI,helixO,bField,xyz[1],dirG,s_out,w_out)) {BadHit(2,tpcHit->position()); continue;}
00420 if (Propagate(lower,normal,helixI,helixO,bField,xyz[2],dirG,s_in ,w_in )) {BadHit(2,tpcHit->position()); continue;}
00421 dx = ((s_out[0] - s_in[0])*w[1] + (s_out[1] - s_in[1])*w[0]);
00422 if (dx <= 0.0) {if (Debug() > 1) {cout << "negative dx " << dx << endl;}
00423 continue;
00424 }
00425
00426 StGlobalDirection globalDirectionOfTrack(dirG);
00427 for (Int_t l = 0; l < 4; l++) {
00428 StGlobalCoordinate globalOfTrack(xyz[l].x(),xyz[l].y(),xyz[l].z());
00429 if (Debug() > 2) {
00430 Int_t k = l;
00431 if (l == 3) k = 0;
00432 Double_t D = - (*mRowPosition[sector-1][row-1][k])*normal;
00433 Double_t A = normal*normal;
00434 Double_t delta = (xyz[l]*normal + D)/TMath::Sqrt(A);
00435 if (TMath::Abs(delta) > 1.e-2) {
00436 cout << "Out of Plane by " << delta << "\tPlane "
00437 << (*mRowPosition[sector-1][row-1][k]) << "\tNormal " << normal << endl;
00438 cout << "Track/hit : " << endl;
00439 cout << "\txyz[0] " << xyz[0] << "\t s = " << s[0] << "/" << s[1] << endl;
00440 cout << "\txyz[1] " << xyz[1] << "\t s_out = " << s_out[0] << "/" << s_out[1] << endl;
00441 cout << "\txyz[2] " << xyz[2] << "\t s_in = " << s_in[0] << "/" << s_in[1] << endl;
00442 cout << "\txyz[3] " << xyz[3] << endl;
00443 }
00444 }
00445 transform(globalOfTrack,localA,sector,row);
00446 transform(localA,localSect[l]);
00447 }
00448 Double_t zP = TMath::Abs(xyz[0].z());
00449
00450 if (zP > 205.0 && zP < 215.) {
00451 Int_t iWestEast = 0;
00452 if (sector > 12) iWestEast = 1;
00453 Int_t io = 0;
00454 if (row > 13) io = 1;
00455 const StThreeVectorD &PromptNormal = *mPromptNormal[iWestEast][io];
00456 const StThreeVectorD &anode = *mPromptPosition[iWestEast][io][0];
00457 const StThreeVectorD &gg = *mPromptPosition[iWestEast][io][1];
00458 const StThreeVectorD &pads = *mPromptPosition[iWestEast][io][2];
00459
00460 if (Propagate(anode,PromptNormal,helixI,helixO,bField,xyz[0],dirG,s,w)) {BadHit(2,tpcHit->position()); continue;}
00461 if (Debug() > 1) {
00462 cout << " Prediction:\t" << xyz[0]
00463 << "\tat s=\t" << s[0] << "/" << s[1]
00464 << "\tw = " << w[0] << "/" << w[1] << endl;
00465 }
00466 dif = xyz[3] - xyz[0];
00467 if (dif.perp() > 2.0) {if (Debug() > 1) {cout << "Prediction is to far from hit:\t" << xyz[3] << endl;}
00468 continue;
00469 }
00470 if (Propagate(pads,PromptNormal,helixI,helixO,bField,xyz[1],dirG,s_outP,w_out)) {BadHit(2,tpcHit->position()); continue;}
00471 if (Propagate(gg ,PromptNormal,helixI,helixO,bField,xyz[2],dirG,s_inP ,w_in )) {BadHit(2,tpcHit->position()); continue;}
00472 s_out[0] = TMath::Min(s_outP[0], s_out[0]);
00473 s_out[1] = TMath::Min(s_outP[1], s_out[1]);
00474 s_in[0] = TMath::Max(s_inP[0] , s_in[0] );
00475 s_in[1] = TMath::Max(s_inP[1] , s_in[1] );
00476 dx = ((s_out[0] - s_in[0])*w[1] + (s_out[1] - s_in[1])*w[0]);
00477 if (dx <= 0.0) {if (Debug() > 1) {cout << "negative dx " << dx << endl;}
00478 continue;
00479 }
00480 StGlobalDirection globalDirectionOfTrack(dirG);
00481 for (Int_t l = 0; l < 4; l++) {
00482 StGlobalCoordinate globalOfTrack(xyz[l].x(),xyz[l].y(),xyz[l].z());
00483 if (Debug() > 2) {
00484 Int_t k = l;
00485 if (l == 3) k = 0;
00486 Double_t D = - (*mRowPosition[sector-1][row-1][k])*PromptNormal;
00487 Double_t A = PromptNormal*PromptNormal;
00488 Double_t delta = (xyz[l]*PromptNormal + D)/TMath::Sqrt(A);
00489 if (TMath::Abs(delta) > 1.e-2) {
00490 cout << "Out of Plane by " << delta << "\tPlane "
00491 << (*mRowPosition[sector-1][row-1][k]) << "\tNormal " << PromptNormal << endl;
00492 cout << "Track/hit : " << endl;
00493 cout << "\txyz[0] " << xyz[0] << "\t s = " << s[0] << "/" << s[1] << endl;
00494 cout << "\txyz[1] " << xyz[1] << "\t s_out = " << s_out[0] << "/" << s_out[1] << endl;
00495 cout << "\txyz[2] " << xyz[2] << "\t s_in = " << s_in[0] << "/" << s_in[1] << endl;
00496 }
00497 }
00498 transform(globalOfTrack,localA,sector,row);
00499 transform(localA,localSect[l]);
00500 }
00501 }
00502 transform(localSect[0],PadOfTrack);
00503 transform(globalDirectionOfTrack,localADirectionOfTrack,sector,row);
00504 transform(localADirectionOfTrack,localDirectionOfTrack);
00505 transform(localSect[3],Pad);
00506 CdEdx[NdEdx].Reset();
00507 CdEdx[NdEdx].resXYZ[0] = localSect[3].position().x() - localSect[0].position().x();
00508 CdEdx[NdEdx].resXYZ[1] = localSect[3].position().y() - localSect[0].position().y();
00509 CdEdx[NdEdx].resXYZ[2] = localSect[3].position().z() - localSect[0].position().z();
00510 TrackLengthTotal += dx;
00511 if (! tpcHit->usedInFit()) {
00512 BadHit(0,tpcHit->position());
00513 continue;}
00514 if ( tpcHit->flag()) {
00515 BadHit(1,tpcHit->position());
00516 continue;}
00517
00518 Int_t iokCheck = 0;
00519 if (sector != Pad.sector() ||
00520 row != Pad.row()) {
00521 gMessMgr->Warning() << "StdEdxY2Maker:: mismatched Sector "
00522 << Pad.sector() << " / " << sector
00523 << " Row " << Pad.row() << " / " << row
00524 << "pad " << Pad.pad() << " TimeBucket :" << Pad.timeBucket()
00525 << endm;
00526 iokCheck++;
00527 }
00528 Double_t pad = tpcHit->pad();
00529 if (pad == 0) pad = Pad.pad();
00530 if (Pad.timeBucket() < 0 ||
00531 Pad.timeBucket() >= numberOfTimeBins) {
00532 gMessMgr->Warning() << "StdEdxY2Maker:: TimeBucket out of range: "
00533 << Pad.timeBucket() << endm;
00534 iokCheck++;
00535 }
00536 if (sector != PadOfTrack.sector() ||
00537 row != PadOfTrack.row() ||
00538 TMath::Abs(Pad.pad()-PadOfTrack.pad()) > 5) {
00539 if (Debug() > 1) {
00540 gMessMgr->Warning() << "StdEdxY2Maker:: Helix Prediction "
00541 << "Sector = "
00542 << PadOfTrack.sector() << "/"
00543 << sector
00544 << " Row = " << PadOfTrack.row() << "/"
00545 << row
00546 << " Pad = " << PadOfTrack.pad() << "/"
00547 << Pad.pad()
00548 << " from Helix is not matched with point/" << endm;;
00549 gMessMgr->Warning() << "StdEdxY2Maker:: Coordinates Preiction: "
00550 << xyz[0] << "/Hit " << tpcHit->position()
00551 << endm;
00552 }
00553 iokCheck++;
00554 }
00555 if (tpcHit->charge() <= 0) {
00556 gMessMgr->Warning() << "StdEdxY2Maker:: deposited charge : " << tpcHit->charge()
00557 << " <= 0" << endm;
00558 iokCheck++;
00559 }
00560
00561 if ((TESTBIT(m_Mode, kPadSelection)) && iokCheck) {BadHit(3, tpcHit->position()); continue;}
00562 if ((TESTBIT(m_Mode, kPadSelection)) && (dx < 0.5 || dx > 25.)) {BadHit(4, tpcHit->position()); continue;}
00563 StTpcdEdxCorrection::ESector kTpcOutIn = StTpcdEdxCorrection::kTpcOuter;
00564 if (row <= NumberOfInnerRows) kTpcOutIn = StTpcdEdxCorrection::kTpcInner;
00565
00566 CdEdx[NdEdx].Reset();
00567 CdEdx[NdEdx].DeltaZ = 5.2;
00568 CdEdx[NdEdx].QRatio = -2;
00569 CdEdx[NdEdx].QRatioA = -2.;
00570 CdEdx[NdEdx].QSumA = 0;
00571 CdEdx[NdEdx].sector = sector;
00572 CdEdx[NdEdx].row = row;
00573 Double_t QAcumm = St_TpcAvgCurrentC::instance()->AcChargeRow(sector,row);
00574 Double_t QL = QAcumm*RA[kTpcOutIn]*RA[kTpcOutIn]/WireLenth[kTpcOutIn];
00575 CdEdx[NdEdx].pad = Pad.pad();
00576 CdEdx[NdEdx].pad = (Int_t) Pad.pad();
00577 CdEdx[NdEdx].edge = CdEdx[NdEdx].pad;
00578 if (CdEdx[NdEdx].edge > 0.5*gStTpcDb->PadPlaneGeometry()->numberOfPadsAtRow(row))
00579 CdEdx[NdEdx].edge += 1 - gStTpcDb->PadPlaneGeometry()->numberOfPadsAtRow(row);
00580 CdEdx[NdEdx].Npads = tpcHit->padsInHit();
00581 CdEdx[NdEdx].Ntbins = tpcHit->pixelsInHit();
00582
00583 CdEdx[NdEdx].dE = tpcHit->charge();
00584
00585 CdEdx[NdEdx].dCharge= tpcHit->chargeModified() - tpcHit->charge();
00586 Int_t p1 = tpcHit->minPad();
00587 Int_t p2 = tpcHit->maxPad();
00588 Int_t t1 = tpcHit->minTmbk();
00589 Int_t t2 = tpcHit->maxTmbk();
00590 CdEdx[NdEdx].rCharge= 0.5*m_TpcdEdxCorrection->Adc2GeV()*TMath::Pi()/4.*(p2-p1+1)*(t2-t1+1);
00591 if (TESTBIT(m_Mode, kEmbeddingShortCut) &&
00592 (tpcHit->idTruth() && tpcHit->qaTruth() > 95)) CdEdx[NdEdx].lSimulated = tpcHit->idTruth();
00593 CdEdx[NdEdx].dx = dx;
00594 CdEdx[NdEdx].dxH = 0;
00595 CdEdx[NdEdx].xyz[0] = localSect[3].position().x();
00596 CdEdx[NdEdx].xyz[1] = localSect[3].position().y();
00597 CdEdx[NdEdx].xyz[2] = localSect[3].position().z();
00598 Double_t probablePad = gStTpcDb->PadPlaneGeometry()->numberOfPadsAtRow(row)/2;
00599 Double_t pitch = (row<14) ?
00600 gStTpcDb->PadPlaneGeometry()->innerSectorPadPitch() :
00601 gStTpcDb->PadPlaneGeometry()->outerSectorPadPitch();
00602 Double_t PhiMax = TMath::ATan2(probablePad*pitch, gStTpcDb->PadPlaneGeometry()->radialDistanceAtRow(row));
00603 CdEdx[NdEdx].PhiR = TMath::ATan2(CdEdx[NdEdx].xyz[0],CdEdx[NdEdx].xyz[1])/PhiMax;
00604 CdEdx[NdEdx].xyzD[0] = localDirectionOfTrack.position().x();
00605 CdEdx[NdEdx].xyzD[1] = localDirectionOfTrack.position().y();
00606 CdEdx[NdEdx].xyzD[2] = localDirectionOfTrack.position().z();
00607 CdEdx[NdEdx].ZdriftDistance = localSect[3].position().z();
00608 CdEdx[NdEdx].zG = tpcHit->position().z();
00609 CdEdx[NdEdx].Qcm = 1e6*QL/tpcHit->position().perp2();
00610 CdEdx[NdEdx].Crow = St_TpcAvgCurrentC::instance()->AvCurrRow(sector,row);
00611 if (St_trigDetSumsC::instance()) CdEdx[NdEdx].Zdc = St_trigDetSumsC::instance()->zdcX();
00612 CdEdx[NdEdx].adc = tpcHit->adc();
00613 Bool_t doIT = kTRUE;
00614 if (TESTBIT(m_Mode,kEmbedding)) doIT = kFALSE;
00615 Int_t iok = m_TpcdEdxCorrection->dEdxCorrection(CdEdx[NdEdx],doIT);
00616 if (iok) {BadHit(4+iok, tpcHit->position()); continue;}
00617 if (fZOfGoodHits) fZOfGoodHits->Fill(tpcHit->position().z());
00618 if (NdEdx < NdEdxMax) {
00619 TrackLength += CdEdx[NdEdx].dx;
00620 NdEdx++;
00621 NoOfTpcHitsUsed++;
00622 }
00623 if (NdEdx > NoFitPoints)
00624 gMessMgr->Error() << "StdEdxY2Maker:: NdEdx = " << NdEdx
00625 << "> NoFitPoints ="<< NoFitPoints << endm;
00626 }
00627 if (fTracklengthInTpcTotal) fTracklengthInTpcTotal->Fill(TrackLengthTotal);
00628 if (fTracklengthInTpc) fTracklengthInTpc->Fill(TrackLength);
00629 SortdEdx();
00630 Double_t I70 = 0, D70 = 0;
00631 Int_t N70 = NdEdx - (int) (0.3*NdEdx + 0.5);
00632 Double_t dXavLog2 = 1;
00633 Int_t k;
00634 #if 0
00635 Double_t SumdEdX = 0;
00636 Double_t SumdX = 0;
00637 for (k = 0; k < NdEdx; k++) {
00638 if (dEdxS[k].dx > 0) {
00639 SumdEdX += dEdxS[k].dEdx;
00640 SumdX += dEdxS[k].dEdx*TMath::Log2(dEdxS[k].dx);
00641 }
00642 }
00643 if (SumdEdX > 0) dXavLog2 = SumdX/SumdEdX;
00644 #endif
00645 if (N70 > 1) {
00646 for (k = 0; k < N70; k++) {
00647 I70 += dEdxS[k].dEdx;
00648 D70 += dEdxS[k].dEdx*dEdxS[k].dEdx;
00649 TrackLength70 += dEdxS[k].dx;
00650 }
00651 I70 /= N70; D70 /= N70;
00652 D70 = TMath::Sqrt(D70 - I70*I70);
00653 D70 /= I70;
00654 dst_dedx_st dedx;
00655 dedx.id_track = Id;
00656 dedx.det_id = kTpcId;
00657 dedx.method = kEnsembleTruncatedMeanId;
00658 dedx.ndedx = N70 + 100*((int) TrackLength);
00659 dedx.dedx[0] = I70;
00660 dedx.dedx[1] = D70;
00661 dedx.dedx[2] = dXavLog2;
00662 if ((TESTBIT(m_Mode, kCalibration)))
00663 for (int l = 0; l < 2; l++) {
00664 if (tracks[l]) tracks[l]->addPidTraits(new StDedxPidTraits(dedx));
00665 }
00666 if (! TESTBIT(m_Mode, kDoNotCorrectdEdx)) {
00667 m_TpcdEdxCorrection->dEdxTrackCorrection(StTpcdEdxCorrection::kTpcLengthCorrection,0,dedx);
00668 m_TpcdEdxCorrection->dEdxTrackCorrection(StTpcdEdxCorrection::kTpcdEdxCor,0,dedx);
00669 dedx.method = kTruncatedMeanId;
00670 for (int l = 0; l < 2; l++) {if (tracks[l]) tracks[l]->addPidTraits(new StDedxPidTraits(dedx));}
00671 }
00672 if ((TESTBIT(m_Mode, kCalibration))) {
00673 Double_t I70A = 0;
00674 Double_t D70A = 0;
00675 Int_t Nouter = 0;
00676 for (k = 0; k < NdEdx; k++) if (dEdxS[k].row > 13) Nouter++;
00677 Int_t N70outer = Nouter - (int) (0.3*Nouter + 0.5);
00678 Double_t TrackLengthA = 0;
00679 if (N70outer > 1) {
00680 Int_t N = 0;
00681 for (k = 0; k < N70outer; k++) {
00682 if (dEdxS[k].row <= 13 || N > N70outer) continue;
00683 N++;
00684 I70A += dEdxS[k].dEdx;
00685 D70A += dEdxS[k].dEdx*dEdxS[k].dEdx;
00686 TrackLengthA += dEdxS[k].dx;
00687 }
00688 if (N > 0) {
00689 I70A /= N; D70A /= N;
00690
00691
00692
00693 D70A = D70A - I70A*I70A;
00694 if (D70A > 0.0){
00695 D70A = TMath::Sqrt(D70A)/I70A;
00696 } else {
00697 D70A = 0;
00698 }
00699 dedx.id_track = Id;
00700 dedx.det_id = kTpcId;
00701 dedx.method = kOtherMethodIdentifier;
00702 dedx.ndedx = N + 100*((int) TrackLengthA);
00703 dedx.dedx[0] = I70A;
00704 dedx.dedx[1] = D70A;
00705 dedx.dedx[2] = dXavLog2;
00706 for (int l = 0; l < 2; l++) {if (tracks[l]) tracks[l]->addPidTraits(new StDedxPidTraits(dedx));}
00707 }
00708 }
00709 }
00710
00711 Double_t chisq, fitZ, fitdZ;
00712 DoFitZ(chisq, fitZ, fitdZ);
00713 if (chisq >0 && chisq < 10000.0) {
00714 dedx.id_track = Id;
00715 dedx.det_id = kTpcId;
00716 dedx.method = kWeightedTruncatedMeanId;
00717 dedx.ndedx = NdEdx + 100*((int) TrackLength);
00718 dedx.dedx[0] = TMath::Exp(fitZ);
00719 dedx.dedx[1] = fitdZ;
00720 dedx.dedx[2] = dXavLog2;
00721 if ((TESTBIT(m_Mode, kCalibration)))
00722 for (int l = 0; l < 2; l++) {if (tracks[l]) tracks[l]->addPidTraits(new StDedxPidTraits(dedx));}
00723 if (! TESTBIT(m_Mode, kDoNotCorrectdEdx)) {
00724 m_TpcdEdxCorrection->dEdxTrackCorrection(StTpcdEdxCorrection::kTpcLengthCorrection,1,dedx);
00725 m_TpcdEdxCorrection->dEdxTrackCorrection(StTpcdEdxCorrection::kTpcdEdxCor,1,dedx);
00726 dedx.method = kLikelihoodFitId;
00727 for (int l = 0; l < 2; l++) {if (tracks[l]) tracks[l]->addPidTraits(new StDedxPidTraits(dedx));}
00728 }
00729 }
00730 if (! TESTBIT(m_Mode, kDoNotCorrectdEdx)) {
00731 StThreeVectorD g3 = gTrack->geometry()->momentum();
00732 Double_t pMomentum = g3.mag();
00733 Float_t Chisq[NHYPS];
00734 for (int hyp = 0; hyp < NHYPS; hyp++) {
00735 Double_t bgL10 = TMath::Log10(pMomentum*TMath::Abs(StProbPidTraits::mPidParticleDefinitions[hyp]->charge())/StProbPidTraits::mPidParticleDefinitions[hyp]->mass());
00736 Chisq[hyp] = LikeliHood(bgL10,NdEdx,FdEdx, StProbPidTraits::mPidParticleDefinitions[hyp]->charge()*StProbPidTraits::mPidParticleDefinitions[hyp]->charge());
00737 }
00738 for (int l = 0; l < 2; l++) {if (tracks[l]) tracks[l]->addPidTraits(new StProbPidTraits(NdEdx,kTpcId,NHYPS,Chisq));}
00739 }
00740 }
00741 }
00742 if ((TESTBIT(m_Mode, kCalibration))) Histogramming(gTrack);
00743 QAPlots(gTrack);
00744 }
00745 if (Debug() > 1) {
00746 gMessMgr->QAInfo() << "StdEdxY2Maker:"
00747 << " Type: " << pEvent->type()
00748 << " Run: " << pEvent->runId()
00749 << " Event: " << pEvent->id()
00750 << " # track nodes: "
00751 << pEvent->trackNodes().size() << endm;
00752 }
00753 if (Debug()) {
00754 timer.stop();
00755 gMessMgr->QAInfo() << "CPU time for StdEdxY2Maker::Make(): "
00756 << timer.elapsedTime() << " sec\n" << endm;
00757 }
00758 if (mHitsUsage) mHitsUsage->Fill(TMath::Log10(TotalNoOfTpcHits+1.), TMath::Log10(NoOfTpcHitsUsed+1.));
00759 return kStOK;
00760 }
00761
00762 void StdEdxY2Maker::SortdEdx() {
00763 Int_t i;
00764 TArrayI idxT(NdEdx); Int_t *idx = idxT.GetArray();
00765 TArrayD dT(NdEdx); Double_t *d = dT.GetArray();
00766 for (i = 0; i < NdEdx; i++) d[i] = CdEdx[i].dEdx;
00767 TMath::Sort(NdEdx,d,idx,0);
00768 for (i=0;i<NdEdx;i++) dEdxS[i] = CdEdx[idx[i]];
00769 TArrayI rowT(NdEdx); Int_t *r = rowT.GetArray();
00770 for (i = 0; i < NdEdx; i++) r[i] = CdEdx[i].row;
00771 TMath::Sort(NdEdx,r,idx,0);
00772 for (i=0;i<NdEdx;i++) FdEdx[i] = CdEdx[idx[i]];
00773 }
00774
00775 void StdEdxY2Maker::Histogramming(StGlobalTrack* gTrack) {
00776 const static Double_t GeV2keV = TMath::Log(1.e-6);
00777
00778 #if 0
00779 static TProfile2D *ETA = 0;
00780 static TH3F *ETA3 = 0;
00781 #endif
00782 static TH2F *Time = 0, *TimeP = 0, *TimeC = 0;
00783 static TH3F *Pressure = 0, *PressureC = 0, *PressureA = 0;
00784 static TH3F *PressureT = 0, *PressureTC = 0, *PressureTA = 0;
00785 static TH3F *Voltage = 0, *VoltageC = 0;
00786 static TH3F *AcCharge = 0;
00787 static TH3F *AvCurrent = 0;
00788 static TH3F *Qcm = 0;
00789
00790 static TH3F **zbgx = 0;
00791
00792 static TH3F *MulRow = 0;
00793 static TH3F *MulRowC = 0;
00794 static TH3F *Phi3 = 0;
00795 static TH3F *Phi3D = 0, *Theta3D = 0;
00796
00797 static TH3F *SecRow3 = 0, *SecRow3C = 0;
00798 static TH3F *Zdc3C = 0;
00799 static TH3F *Z3C = 0, *Z3OC = 0;
00800 #if 0
00801 static TH3F *Z3 = 0, *Z3A = 0;
00802 static TH3F *Z3O = 0, *Z3OA = 0;
00803 static TH3F *Z3OW = 0, *Z3OWC = 0, *Z3OWA = 0;
00804 #endif
00805
00806 static TH2F *AdcI = 0, *AdcO = 0, *AdcIC = 0, *AdcOC = 0, *Adc3I = 0, *Adc3O = 0, *Adc3IC = 0, *Adc3OC = 0;
00807 static TH3F *AdcIZP = 0, *AdcOZP = 0, *AdcIZN = 0, *AdcOZN = 0;
00808 static TH2F **Adc3Ip = 0, **Adc3Op = 0;
00809
00810 static TH2F *ZdcCP = 0, *BBCP = 0, *MultiplicityPI = 0, *MultiplicityPO = 0;
00811 static TH3F *BBC3 = 0;
00812
00813 static TH2F *bbcYellowBkg = 0, *bbcBlueBkg = 0;
00814
00815 static TH3F *SecRow3Mip = 0;
00816
00817
00818 static TH1F *hdEI = 0, *hdEUI = 0, *hdERI = 0, *hdEPI = 0, *hdETI = 0, *hdESI = 0, *hdEZI = 0, *hdEMI = 0;
00819 static TH1F *hdEO = 0, *hdEUO = 0, *hdERO = 0, *hdEPO = 0, *hdETO = 0, *hdESO = 0, *hdEZO = 0, *hdEMO = 0;
00820 #if 0
00821
00822
00823
00824 static TH2F *Points[25][5];
00825 #endif
00826 static TH2F *TPoints[25][5];
00827 #if 0
00828 static TH2F *MPoints[25][5];
00829 #endif
00830 static TH2F *hist70B[NHYPS][2], *histzB[NHYPS][2];
00831 static TH2F *hist70BT[NHYPS][2], *histzBT[NHYPS][2];
00832 static TProfile *hitsB[NHYPS][2];
00833 static TH2F *FitPull = 0, *Pull70 = 0;
00834 static TTree *ftree = 0;
00835 static TH2F *ffitZ[NHYPS], *ffitP[NHYPS], *ffitZU = 0, *ffitZU3 = 0, *ffitZA = 0;
00836 const static Int_t Nlog2dx = 140;
00837 const static Double_t log2dxLow = 0.0, log2dxHigh = 3.5;
00838 static TH2F *inputTPCGasPressureP = 0, *nitrogenPressureP = 0, *gasPressureDiffP = 0, *inputGasTemperatureP = 0;
00839 static TH2F *outputGasTemperatureP = 0, *flowRateArgon1P = 0, *flowRateArgon2P = 0;
00840 static TH2F *flowRateMethaneP = 0;
00841 static TH2F *percentMethaneInP = 0, *percentMethaneInPC = 0, *percentMethaneInPA = 0;
00842 static TH2F *ppmOxygenInP = 0, *flowRateExhaustP = 0;
00843 static TH2F *flowRateRecirculationP = 0;
00844 static TH2F *ppmWaterOutP = 0, *ppmWaterOutPC = 0, *ppmWaterOutPA = 0;
00845
00846 static TH3F *Prob = 0;
00847
00848 static TH3F *dXdE = 0, *dXdEA = 0, *dXdEC = 0;
00849
00850 static TH2F *BaddEdxZPhi70[2], *BaddEdxZPhiZ[2];
00851 static TH1F *BaddEdxMult70[2], *BaddEdxMultZ[2];
00852 static TH3D *Edge3 = 0;
00853 static dEdxTrackY2 *ftrack = 0;
00854 static int hMade = 0;
00855
00856 if (! gTrack && !hMade) {
00857 hMade=2004;
00858
00859 Int_t nZBins = 200;
00860 Double_t ZdEdxMin = -5;
00861 Double_t ZdEdxMax = 5;
00862 #if 0
00863 Z3 = new TH3F("Z3",
00864 "log(dEdx/Pion) versus row and Drift Distance",
00865 NumberOfRows,0.5, NumberOfRows+0.5,105,0.,210.,nZBins,ZdEdxMin,ZdEdxMax);
00866 Z3A = new TH3F("Z3A",
00867 "log(dEdx/Pion) just after corection versus row and Drift Distance",
00868 NumberOfRows,0.5, NumberOfRows+0.5,105,0.,210.,nZBins,ZdEdxMin,ZdEdxMax);
00869 #endif
00870 Z3C = new TH3F("Z3C",
00871 "log(dEdx/Pion) corrected versus row and Drift Distance",
00872 NumberOfRows,0.5, NumberOfRows+0.5,105,0.,210.,nZBins,ZdEdxMin,ZdEdxMax);
00873 #if 0
00874 Z3O = new TH3F("Z3O",
00875 "log(dEdx/Pion) versus row and (Drift)*ppmO2In",
00876 NumberOfRows,0.5, NumberOfRows+0.5,100,0.,1.e4,nZBins,ZdEdxMin,ZdEdxMax);
00877 Z3OA= new TH3F("Z3OA",
00878 "log(dEdx/Pion) just after correction versus row and (Drift)*ppmO2In",
00879 NumberOfRows,0.5, NumberOfRows+0.5,100,0.,1.e4,nZBins,ZdEdxMin,ZdEdxMax);
00880 #endif
00881 Z3OC= new TH3F("Z3OC",
00882 "log(dEdx/Pion) corrected versus row and (Drift)*ppmO2In",
00883 NumberOfRows,0.5, NumberOfRows+0.5,100,0.,1.e4,nZBins,ZdEdxMin,ZdEdxMax);
00884 #if 0
00885 Z3OW = new TH3F("Z3OW",
00886 "log(dEdx/Pion) versus row and (Drift)*ppmO2In*ppmWaterOut*",
00887 NumberOfRows,0.5, NumberOfRows+0.5,100,0.,1.2e5,nZBins,ZdEdxMin,ZdEdxMax);
00888 Z3OWA= new TH3F("Z3OWA",
00889 "log(dEdx/Pion) just after correction versus row and (Drift)*ppmO2In*ppmWaterOut",
00890 NumberOfRows,0.5, NumberOfRows+0.5,100,0.,1.2e5,nZBins,ZdEdxMin,ZdEdxMax);
00891 Z3OWC= new TH3F("Z3OWC",
00892 "log(dEdx/Pion) corrected versus row and (Drift)*ppmO2In*ppmWaterOut",
00893 NumberOfRows,0.5, NumberOfRows+0.5,100,0.,1.2e5,nZBins,ZdEdxMin,ZdEdxMax);
00894 #endif
00895 Edge3 = new TH3D("Edge3",
00896 "log(dEdx/Pion) versus row and Edge",
00897 NumberOfRows,0.5, NumberOfRows+0.5, 400,-100,100,nZBins,ZdEdxMin,ZdEdxMax);
00898 #if 0
00899
00900 ETA = new TProfile2D("ETA",
00901 "log(dEdx/Pion) versus Sector I/O and #{eta}",
00902 NumberOfRows,0.5, NumberOfRows+0.5, 135,-2.25,2.25);
00903 ETA3 = new TH3F("ETA3",
00904 "log(dEdx/Pion) versus Sector I/O and #{eta}",
00905 NumberOfRows,0.5, NumberOfRows+0.5, 135,-2.25,2.25,nZBins,ZdEdxMin,ZdEdxMax);
00906 #endif
00907 SecRow3= new TH3F("SecRow3","<log(dEdx/Pion)> (uncorrected) versus sector and row",
00908 numberOfSectors,0.5, numberOfSectors+0.5, NumberOfRows,0.5, NumberOfRows+0.5,nZBins,ZdEdxMin,ZdEdxMax);
00909 SecRow3->SetXTitle("Sector number");
00910 SecRow3->SetYTitle("Row number");
00911 SecRow3C= new TH3F("SecRow3C","<log(dEdx/Pion)> (corrected) versus sector and row",
00912 numberOfSectors,0.5, numberOfSectors+0.5, NumberOfRows,0.5, NumberOfRows+0.5,nZBins,ZdEdxMin,ZdEdxMax);
00913 SecRow3C->SetXTitle("Sector number");
00914 SecRow3C->SetYTitle("Row number");
00915 #if 0
00916 SecRow3A= new TH3F("SecRow3A","<log(dEdx/Pion)> (just after correction) versus sector and row",
00917 numberOfSectors,0.5, numberOfSectors+0.5, NumberOfRows,0.5, NumberOfRows+0.5,nZBins,ZdEdxMin,ZdEdxMax);
00918 SecRow3A->SetXTitle("Sector number");
00919 SecRow3A->SetYTitle("Row number");
00920 #endif
00921 Zdc3C = new TH3F("Zdc3C","<log(dEdx/Pion)> versus row and ZdcCoincidenceRate (log10)",
00922 NumberOfRows,0.5, NumberOfRows+0.5,100,0,10,nZBins,ZdEdxMin,ZdEdxMax);
00923 MultiplicityPI = new TH2F("MultiplicityPI","Multiplicity (log10) Inner",100,0,10,nZBins,ZdEdxMin,ZdEdxMax);
00924 MultiplicityPO = new TH2F("MultiplicityPO","Multiplicity (log10) Outer",100,0,10,nZBins,ZdEdxMin,ZdEdxMax);
00925 if ((TESTBIT(m_Mode, kAdcHistos))) {
00926 gMessMgr->Warning() << "StdEdxY2Maker::Histogramming Adc Histograms" << endm;
00927 AdcI = new TH2F("AdcI",
00928 "log10dE (keV measured) versus log10dE(Predicted) for Inner rows",
00929 120,-5.7,-3.3,185,-7.2,-3.5);
00930 AdcO = new TH2F("AdcO",
00931 "log10dE (keV measured) versus log10dE(Predicted) for Outer rows",
00932 150,-5.5,-2.5,165,-6.5,-3.0);
00933 AdcIC = new TH2F("AdcIC",
00934 "log10dE (keV measured corrected) versus log10dE(Predicted) for Inner rows",
00935 120,-5.7,-3.3,185,-7.2,-3.5);
00936 AdcOC = new TH2F("AdcOC",
00937 "log10dE (keV measured corrected) versus log10dE(Predicted) for Outer rows",
00938 150,-5.5,-2.5,165,-6.5,-3.0);
00939 Adc3I = new TH2F("Adc3I",
00940 "Uniq 3*sigma log10dE (keV measured) versus log10dE(Predicted) for Inner rows",
00941 120,-5.7,-3.3,185,-7.2,-3.5);
00942 Adc3O = new TH2F("Adc3O",
00943 "Uniq 3*sigma log10dE (keV measured) versus log10dE(Predicted) for Outer rows",
00944 150,-5.5,-2.5,165,-6.5,-3.0);
00945 Adc3IC = new TH2F("Adc3IC",
00946 "Uniq 3*sigma log10dE (keV measured corrected) versus log10dE(Predicted) for Inner rows",
00947 120,-5.7,-3.3,185,-7.2,-3.5);
00948 Adc3OC = new TH2F("Adc3OC",
00949 "Uniq 3*sigma log10dE (keV measured corrected) versus log10dE(Predicted) for Outer rows",
00950 150,-5.5,-2.5,165,-6.5,-3.0);
00951 Adc3Ip = new TH2F*[NHYPS];
00952 Adc3Op = new TH2F*[NHYPS];
00953 for (int hyp = 0; hyp < NHYPS; hyp++) {
00954 TString nameP(StProbPidTraits::mPidParticleDefinitions[hyp]->name().data());
00955 nameP.ReplaceAll("-","");
00956 Adc3Ip[hyp] = new
00957 TH2F(Form("Adc3I%s",nameP.Data()),
00958 Form("%s Uniq 3*sigma log10dE (keV meas.cor.) versus log10dE(Predicted) for Inner rows",
00959 nameP.Data()),
00960 120,-5.7,-3.3,185,-7.2,-3.5);
00961 Adc3Op[hyp] = new
00962 TH2F(Form("Adc3O%s",nameP.Data()),
00963 Form("%s Uniq 3*sigma log10dE (keV meas.cor.) versus log10dE(Predicted) for Outer rows",
00964 nameP.Data()),
00965 120,-5.7,-3.3,185,-7.2,-3.5);
00966 }
00967 AdcIZP = new TH3F("AdcIZP","z (Positive measured) versus dE(Predicted) and Z for Inner rows",
00968 200,0,200,100,0.,1.e4,nZBins,ZdEdxMin,ZdEdxMax);
00969 AdcOZP = new TH3F("AdcOZP","z (Positive measured) versus dE(Predicted) and Z for Outer rows",
00970 500,0,500,100,0.,1.e4,nZBins,ZdEdxMin,ZdEdxMax);
00971 AdcIZN = new TH3F("AdcIZN","z (Positive measured) versus dE(Predicted) and Z for Inner rows",
00972 200,0,200,100,0.,1.e4,nZBins,ZdEdxMin,ZdEdxMax);
00973 AdcOZN = new TH3F("AdcOZN","z (Positive measured) versus dE(Predicted) and Z for Outer rows",
00974 500,0,500,100,0.,1.e4,nZBins,ZdEdxMin,ZdEdxMax);
00975 }
00976 ZdcCP = new TH2F("ZdcCP","ZdcCoincidenceRate (log10)",100,0,10,nZBins,ZdEdxMin,ZdEdxMax);
00977 BBCP = new TH2F("BBCP","BbcCoincidenceRate (log10)",60,0,6,nZBins,ZdEdxMin,ZdEdxMax);
00978 BBC3 = new TH3F("BBC3","BbcCoincidenceRate (log10) and row ",NumberOfRows,0.5, NumberOfRows+0.5,60,0,6,nZBins,ZdEdxMin,ZdEdxMax);
00979 bbcYellowBkg = new TH2F("bbcYellowBkg","(BBC Eastdelayed) and (BBC West) (log10)",
00980 100,0,10,nZBins,ZdEdxMin,ZdEdxMax);
00981 bbcBlueBkg = new TH2F("bbcBlueBkg","(BBC Westdelayed) and (BBC East) (log10)",
00982 100,0,10,nZBins,ZdEdxMin,ZdEdxMax);
00983 if ((TESTBIT(m_Mode, kMip))) {
00984 gMessMgr->Warning() << "StdEdxY2Maker::Histogramming Mip Histograms" << endm;
00985 SecRow3Mip = new TH3F
00986 ("SecRow3Mip",
00987 "<log(dEdx/Pion)>/sigma (corrected) versus row and log2(dx) for MIP particle)",
00988 NumberOfRows,0.5, NumberOfRows+0.5,Nlog2dx, log2dxLow, log2dxHigh, 200,-5,15);
00989 }
00990 MulRow = new TH3F("MulRow","log(dEdx/Pion) versus log10(Multplicity) and row",
00991 100,0,10, NumberOfRows,0.5, NumberOfRows+0.5,nZBins,ZdEdxMin,ZdEdxMax);
00992 MulRowC = new TH3F("MulRowC","log(dEdx/Pion) versus log10(Multplicity) and row corrected",
00993 100,0,10, NumberOfRows,0.5, NumberOfRows+0.5,nZBins,ZdEdxMin,ZdEdxMax);
00994 Phi3 = new TH3F("Phi3","log(dEdx/Pion) versus Phi (coordinates, relative) (degrees) and row",
00995 210,-1.05,1.05, NumberOfRows,0.5, NumberOfRows+0.5,nZBins,ZdEdxMin,ZdEdxMax);
00996 Phi3D = new TH3F("Phi3D","log(dEdx/Pion) versus Phi (direction) and row",
00997 480,-60,60, NumberOfRows,0.5, NumberOfRows+0.5,nZBins,ZdEdxMin,ZdEdxMax);
00998 Theta3D = new TH3F("Theta3D","log(dEdx/Pion) versus Theta (direction) (degrees) and row",
00999 560,-60,80, NumberOfRows,0.5, NumberOfRows+0.5,nZBins,ZdEdxMin,ZdEdxMax);
01000 const Char_t *N[5] = {"B","70B","BU","70BU","70BA"};
01001 const Char_t *T[5] = {"dEdx(fit)/Pion",
01002 "dEdx(fit_uncorrected)/Pion ",
01003 "dEdx(I70)/Pion",
01004 "dEdx(I70_uncorrected)/Pion",
01005 "dEdx(I70A_uncorrected)/Pion"};
01006
01007 Int_t NZ = 1;
01008 for (Int_t t = 0; t < 5; t++) {
01009 for (Int_t z = 0; z < NZ; z++) {
01010 TString ZN("");
01011 TString ZT("all");
01012 if (z > 0) {
01013 ZN = Form("%02i",z);
01014 ZT = Form("Sector %02i",z);
01015 }
01016 #if 0
01017 Points[z][t] = new TH2F(Form("Points%s%s",N[t],ZN.Data()),
01018 Form("%s/sigma versus no. of measured points for %s",T[t],ZT.Data()),
01019 50,0,50., 500,-5.,20.);
01020 MPoints[z][t] = new TH2F(Form("MPoints%s%s",N[t],ZN.Data()),
01021 Form("%s versus no. of measured points for %s",T[t],ZT.Data()),
01022 150,10,160., 500,-1.,4.);
01023 #endif
01024 TPoints[z][t] = new TH2F(Form("TPoints%s%s",N[t],ZN.Data()),
01025 Form("%s versus Length in Tpc for %s",T[t],ZT.Data()),
01026 150,10,160., 500,-1.,4.);
01027 }
01028 }
01029 if ((TESTBIT(m_Mode, kZBGX))) {
01030 gMessMgr->Warning() << "StdEdxY2Maker::Histogramming make zbgx histograms" << endm;
01031 zbgx = new TH3F*[2*NHYPS];
01032 }
01033 for (int hyp=0; hyp<NHYPS;hyp++) {
01034 TString nameP("fit");
01035 nameP += StProbPidTraits::mPidParticleDefinitions[hyp]->name().data();
01036 nameP.ReplaceAll("-","");
01037 TString title = "fitZ - Pred. for ";
01038 title += StProbPidTraits::mPidParticleDefinitions[hyp]->name().data();
01039 title.ReplaceAll("-","");
01040 title += " versus log10(beta*gamma) for pion";
01041 ffitZ[hyp] = new TH2F(nameP.Data(),title.Data(),100,-1,4,100,-5,5);
01042 ffitZ[hyp]->SetMarkerColor(hyp+2);
01043 ffitP[hyp] = new TH2F(*ffitZ[hyp]);
01044 nameP.ReplaceAll("fit","fitP");
01045 ffitP[hyp]->SetName(nameP.Data());
01046 for (int sCharge = 0; sCharge < 2; sCharge++) {
01047 if ((TESTBIT(m_Mode, kZBGX))) {
01048 nameP = "zbgx";
01049 nameP += StProbPidTraits::mPidParticleDefinitions[hyp]->name().data();
01050 nameP.ReplaceAll("-","");
01051 if (sCharge == 0) nameP += "P";
01052 else nameP += "N";
01053 if (zbgx)
01054 zbgx[2*hyp+sCharge] = new TH3F(nameP.Data(),"z = log(dE/dx) versus log10(beta*gamma) and log2(dx) for unique hyps",
01055 100,-1,4,Nlog2dx,log2dxLow,log2dxHigh,320,-2,6);
01056 }
01057 nameP = StProbPidTraits::mPidParticleDefinitions[hyp]->name().data();
01058 nameP.ReplaceAll("-","");
01059 if (sCharge == 0) nameP += "P";
01060 else nameP += "N";
01061 TString name = nameP;
01062 name += "70";
01063 title = "log(dE/dx70/I(";
01064 title += nameP;
01065 title += ")) versus log10(p/m)";
01066 name += "B";
01067 title += " Bichsel";
01068 hist70B[hyp][sCharge] = new TH2F(name.Data(),title.Data(),100,-1.,4.,600,-2.,4.);
01069 name += "T";
01070 title += " Unique";
01071 hist70BT[hyp][sCharge] = new TH2F(name.Data(),title.Data(),100,-1.,4.,600,-2.,4.);
01072 name = nameP;
01073 name += "z";
01074 title = "zFit - log(I(";
01075 title += nameP;
01076 title += ")) versus log10(p/m)";
01077 name += "B";
01078 title += " Bichsel";
01079 histzB[hyp][sCharge] = new TH2F(name.Data(),title.Data(),100,-1.,4.,600,-2.,4.);
01080 name += "T";
01081 title += " Unique";
01082 histzBT[hyp][sCharge] = new TH2F(name.Data(),title.Data(),100,-1.,4.,600,-2.,4.);
01083 name = nameP;
01084 name += "B";
01085 name += "B";
01086 title = "log(I_{BB}(";
01087 title += nameP;
01088 title += ")) versus log10(p/m) Bichsel";
01089 hitsB[hyp][sCharge] = new TProfile(name.Data(),title.Data(),100,-1.,4.);
01090 }
01091 }
01092 TDatime t1(tMin,0);
01093 TDatime t2(tMax,0);
01094
01095 UInt_t i1 = t1.Convert() - timeOffSet;
01096 UInt_t i2 = t2.Convert() - timeOffSet;
01097 Int_t Nt = (i2 - i1)/(3600);
01098 Pressure = new TH3F("Pressure","log(dE/dx)_{uncorrected} - log(I(pi)) versus Row & Log(Pressure)",
01099 NumberOfRows,0.5, NumberOfRows+0.5,150, 6.84, 6.99,nZBins,ZdEdxMin,ZdEdxMax);
01100 PressureA = new TH3F("PressureA","log(dE/dx)_{just after correction} log(I(pi)) versus Log(Pressure)",
01101 NumberOfRows,0.5, NumberOfRows+0.5,150, 6.84, 6.99,nZBins,ZdEdxMin,ZdEdxMax);
01102 PressureC = new TH3F("PressureC","log(dE/dx)_{corrected} - row & log(I(pi)) versus Log(Pressure)",
01103 NumberOfRows,0.5, NumberOfRows+0.5,150, 6.84, 6.99,nZBins,ZdEdxMin,ZdEdxMax);
01104 PressureT = new TH3F("PressureT","log(dE/dx)_{uncorrected} - log(I(pi)) versus Row & Log(Pressure*298.2/outputGasTemperature)",
01105 NumberOfRows,0.5, NumberOfRows+0.5,150, 6.84, 6.99,nZBins,ZdEdxMin,ZdEdxMax);
01106 PressureTA = new TH3F("PressureTA","log(dE/dx)_{just after correction} log(I(pi)) versus Log(Pressure*298.2/outputGasTemperature)",
01107 NumberOfRows,0.5, NumberOfRows+0.5,150, 6.84, 6.99,nZBins,ZdEdxMin,ZdEdxMax);
01108 PressureTC = new TH3F("PressureTC","log(dE/dx)_{corrected} - row & log(I(pi)) versus Log(Pressure*298.2/outputGasTemperature)",
01109 NumberOfRows,0.5, NumberOfRows+0.5,150, 6.84, 6.99,nZBins,ZdEdxMin,ZdEdxMax);
01110 Voltage = new TH3F("Voltage","log(dE/dx)_{uncorrected} - log(I(pi)) versus Sector*Row and Voltage - Voltage_{nominal}",
01111 numberOfSectors*NumberOfRows,0.5, numberOfSectors*NumberOfRows+0.5,6,-110,10,nZBins,ZdEdxMin,ZdEdxMax);
01112 VoltageC = new TH3F("VoltageC","log(dE/dx)_{corrected} - row & log(I(pi)) versus Row and Voltage",
01113 NumberOfRows,0.5, NumberOfRows+0.5,410,990.,1400.,nZBins,ZdEdxMin,ZdEdxMax);
01114 AcCharge = new TH3F("AcCharge","log(dE/dx)_{corrected} - row & log(I(pi)) versus Row and Accumulated Charge [C]",
01115 NumberOfRows,0.5, NumberOfRows+0.5,110,0.,5.,nZBins,ZdEdxMin,ZdEdxMax);
01116 AvCurrent = new TH3F("AvCurrent","log(dE/dx)_{corrected} - row & log(I(pi)) versus Row and Average Current [#{mu}A]",
01117 NumberOfRows,0.5, NumberOfRows+0.5,210,0.0,1.0,nZBins,ZdEdxMin,ZdEdxMax);
01118 Qcm = new TH3F("Qcm","log(dE/dx)_{corrected} - row & log(I(pi)) versus Row and Accumulated Charge [uC/cm]",
01119 NumberOfRows,0.5, NumberOfRows+0.5,100,0.,10.,nZBins,ZdEdxMin,ZdEdxMax);
01120
01121
01122 Time = new TH2F("Time","log(dE/dx)_{uncorrected} - log(I(pi)) versus Date& Time",
01123 Nt,i1,i2,nZBins,ZdEdxMin,ZdEdxMax);
01124 TimeC = new TH2F("TimeC","log(dE/dx)_{corrected} - log(I(pi)) versus Date& Time after correction",
01125 Nt,i1,i2,nZBins,ZdEdxMin,ZdEdxMax);
01126 TimeP = new TH2F("TimeP","log(dE/dx)_{after pressure correction} - log(I(pi)) versus Date& Time",
01127 Nt,i1,i2,nZBins,ZdEdxMin,ZdEdxMax);
01128 FitPull= new TH2F("FitPull","(zFit - log(I(pi)))/dzFit versus track length",
01129 150,10.,160,nZBins,ZdEdxMin,ZdEdxMax);
01130 Pull70 = new TH2F("Pull70","log(I70/I(pi)))/D70 versus track length",
01131 150,10.,160,nZBins,ZdEdxMin,ZdEdxMax);
01132 TString title("");
01133 title = "log(dE/dx/Pion) vs inputTPCGasPressure (mbar)";
01134 inputTPCGasPressureP = new TH2F("inputTPCGasPressureP","log(dE/dx/Pion) vs inputTPCGasPressure (mbar)",100,2.0,2.2,nZBins,ZdEdxMin,ZdEdxMax);
01135 nitrogenPressureP = new TH2F("nitrogenPressureP","log(dE/dx/Pion) vs nitrogenPressure (mbar)",100,0.9,1.1,nZBins,ZdEdxMin,ZdEdxMax);
01136 gasPressureDiffP = new TH2F("gasPressureDiffP","log(dE/dx/Pion) vs gasPressureDiff (mbar)",100,0.6,1.,nZBins,ZdEdxMin,ZdEdxMax);
01137 inputGasTemperatureP = new TH2F("inputGasTemperatureP","log(dE/dx/Pion) vs inputGasTemperature (degrees C)",100,295.,300.,nZBins,ZdEdxMin,ZdEdxMax);
01138 outputGasTemperatureP = new TH2F("outputGasTemperatureP","log(dE/dx/Pion) vs outputGasTemperature (degrees C)",100,295.,300.,nZBins,ZdEdxMin,ZdEdxMax);
01139 flowRateArgon1P = new TH2F("flowRateArgon1P","log(dE/dx/Pion) vs flowRateArgon1 (liters/min)",100,14.95,15.0,nZBins,ZdEdxMin,ZdEdxMax);
01140 flowRateArgon2P = new TH2F("flowRateArgon2P","log(dE/dx/Pion) vs flowRateArgon2 (liters/min)",100,0.,0.25,nZBins,ZdEdxMin,ZdEdxMax);
01141 flowRateMethaneP = new TH2F("flowRateMethaneP","log(dE/dx/Pion) vs flowRateMethane (liters/min)",100,1.34,1.37,nZBins,ZdEdxMin,ZdEdxMax);
01142 percentMethaneInP = new TH2F("percentMethaneInP","log(dE/dx/Pion) vs percentMethaneIn (percent)",100,9.6,10.6,nZBins,ZdEdxMin,ZdEdxMax);
01143 percentMethaneInPC = new TH2F("percentMethaneInPC","log(dE/dx/Pion)(corrected) vs percentMethaneIn (percent)",100,9.6,10.6,nZBins,ZdEdxMin,ZdEdxMax);
01144 percentMethaneInPA = new TH2F("percentMethaneInPA","log(dE/dx/Pion)(just after correction) vs percentMethaneIn (percent)",
01145 100,9.6,10.6,nZBins,ZdEdxMin,ZdEdxMax);
01146 ppmOxygenInP = new TH2F("ppmOxygenInP","log(dE/dx/Pion) vs ppmOxygenIn (ppm)",240,0.,60.,nZBins,ZdEdxMin,ZdEdxMax);
01147 flowRateExhaustP = new TH2F("flowRateExhaustP","log(dE/dx/Pion) vs flowRateExhaust (liters/min)",100,5.,20.,nZBins,ZdEdxMin,ZdEdxMax);
01148 ppmWaterOutP = new TH2F("ppmWaterOutP","log(dE/dx/Pion) vs ppmWaterOut (ppm)",100,0.,20.,nZBins,ZdEdxMin,ZdEdxMax);
01149 ppmWaterOutPC = new TH2F("ppmWaterOutPC","log(dE/dx/Pion) corrected vs ppmWaterOut (ppm)",100,0.,20.,nZBins,ZdEdxMin,ZdEdxMax);
01150 ppmWaterOutPA = new TH2F("ppmWaterOutPA","log(dE/dx/Pion) just after correction vs ppmWaterOut (ppm)",100,0.,20.,nZBins,ZdEdxMin,ZdEdxMax);
01151
01152 flowRateRecirculationP = new TH2F("flowRateRecirculationP","log(dE/dx/Pion) vs flowRateRecirculation (liters/min)",
01153 100,515.,545.,nZBins,ZdEdxMin,ZdEdxMax);
01154 ffitZU = new TH2F("fitZU","fitZ - PredPi Unique versus log10(beta*gamma)",100,-1,4,100,-5,5);
01155 ffitZU->SetMarkerColor(7);
01156 ffitZU3 = new TH2F("fitZU3","fitZ - PredPi Unique and 3 sigma away versus log10(beta*gamma)",100,-1,4,100,-5,5);
01157 ffitZU3->SetMarkerColor(6);
01158 ffitZA = new TH2F("fitZA","fitZ - PredPi All versus log10(beta*gamma)",100,-1,4,100,-5,5);
01159 ffitZA->SetMarkerColor(1);
01160 hdEI = new TH1F("hdEI","log10(dE) Inner after calibration",100,-8.,-3.);
01161 hdEUI = new TH1F("hdEUI","log10(dEU) Inner before correction",100,-8.,-3.);
01162 hdERI = new TH1F("hdERI","log10(dER) Inner after row correction correction",100,-8.,-3.);
01163 hdEPI = new TH1F("hdEPI","log10(dEP) Inner after Pressure correction",100,-8.,-3.);
01164 hdETI = new TH1F("hdETI","log10(dET) Inner after TimeScale",100,-8.,-3.);
01165 hdESI = new TH1F("hdESI","log10(dES) Inner after after TimeScale + SecRow corrections",100,-8.,-3.);
01166 hdEZI = new TH1F("hdEZI","log10(dEZ) Inner after TimeScale + SecRow + Sec Z corrections ",100,-8.,-3.);
01167 hdEMI = new TH1F("hdEMI","log10(dEM) Inner after TimeScale + SecRow + Sec Z + Multiplicity corrections",100,-8.,-3.);
01168 hdEO = new TH1F("hdEO","log10(dE) Outer after calibration",100,-8.,-3.);
01169 hdEUO = new TH1F("hdEUO","log10(dEU) Outer before correction",100,-8.,-3.);
01170 hdERO = new TH1F("hdERO","log10(dER) Outer after row correction correction",100,-8.,-3.);
01171 hdEPO = new TH1F("hdEPO","log10(dEP) Outer after Pressure correction",100,-8.,-3.);
01172 hdETO = new TH1F("hdETO","log10(dET) Outer after TimeScale",100,-8.,-3.);
01173 hdESO = new TH1F("hdESO","log10(dES) Outer after after TimeScale + SecRow corrections",100,-8.,-3.);
01174 hdEZO = new TH1F("hdEZO","log10(dEZ) Outer after TimeScale + SecRow + Sec Z corrections ",100,-8.,-3.);
01175 hdEMO = new TH1F("hdEMO","log10(dEM) Outer after TimeScale + SecRow + Sec Z + Multiplicity corrections",100,-8.,-3.);
01176 if ((TESTBIT(m_Mode, kProbabilityPlot))) {
01177 gMessMgr->Warning() << "StdEdxY2Maker::Histogramming Probability Histograms" << endm;
01178 Prob = new TH3F("Prob","Z(=log(I70/Bichsel)) versun log10(bg) for pion and Probability",
01179 100,-1.,4.,10*NHYPS+1,-.1,NHYPS,600,-2.,4.);
01180 }
01181 dXdE = new TH3F("dXdE","log(dEdx/Pion) versus dX and row",
01182 100,0.,5., NumberOfRows,0.5, NumberOfRows+0.5,nZBins,ZdEdxMin,ZdEdxMax);
01183 dXdEA = new TH3F("dXdEA","log(dEdx/Pion) just after correction versus dX and row",
01184 100,0.,5., NumberOfRows,0.5, NumberOfRows+0.5,nZBins,ZdEdxMin,ZdEdxMax);
01185 dXdEC = new TH3F("dXdEC","log(dEdx/Pion) corrected versus dX and row",
01186 100,0.,5., NumberOfRows,0.5, NumberOfRows+0.5,nZBins,ZdEdxMin,ZdEdxMax);
01187
01188
01189 BaddEdxZPhi70[0] = new TH2F("BaddEdxZPhi700","Z and Phi for I70 below any limits by 5 s.d.",210,-210,210,360,-180.,180.);
01190 BaddEdxZPhi70[1] = new TH2F("BaddEdxZPhi701","Z and Phi for I70 above any limits by 5 s.d.",210,-210,210,360,-180.,180.);
01191 BaddEdxMult70[0] = new TH1F("BaddEdxMult700","Multiplicity (log10) for I70 below any limits by 5 s.d.",100,0.,10.);
01192 BaddEdxMult70[1] = new TH1F("BaddEdxMult701","Multiplicity (log10) for I70 above any limits by 5 s.d.",100,0.,10.);
01193 BaddEdxZPhiZ[0] = new TH2F("BaddEdxZPhiZ0","Z and Phi for Ifit below any limits by 5 s.d.",210,-210,210,360,-180.,180.);
01194 BaddEdxZPhiZ[1] = new TH2F("BaddEdxZPhiZ1","Z and Phi for Ifit above any limits by 5 s.d.",210,-210,210,360,-180.,180.);
01195 BaddEdxMultZ[0] = new TH1F("BaddEdxMultZ0","Multiplicity (log10) for Ifit below any limits by 5 s.d.",100,0.,10.);
01196 BaddEdxMultZ[1] = new TH1F("BaddEdxMultZ1","Multiplicity (log10) for Ifit above any limits by 5 s.d.",100,0.,10.);
01197 if ((TESTBIT(m_Mode, kMakeTree))&& !ftree) {
01198 gMessMgr->Warning() << "StdEdxY2Maker::Histogramming Make Tree" << endm;
01199 ftree = new TTree("dEdxT","dEdx tree");
01200 ftree->SetAutoSave(1000000000);
01201 Int_t bufsize = 64000;
01202 Int_t split = 99;
01203 if (split) bufsize /= 4;
01204 ftrack = new dEdxTrackY2();
01205 TTree::SetBranchStyle(1);
01206 TBranch *branch = ftree->Branch("dEdxTrackY2", "dEdxTrackY2", &ftrack, bufsize,split);
01207 branch->SetAutoDelete(kFALSE);
01208 }
01209 return;
01210 }
01211
01212 tpcGas_st *tpcGas = 0;
01213 if ( m_TpcdEdxCorrection && m_TpcdEdxCorrection->tpcGas()) tpcGas = m_TpcdEdxCorrection->tpcGas()->GetTable();
01214
01215 StThreeVectorD g3 = gTrack->geometry()->momentum();
01216 Double_t pMomentum = g3.mag();
01217 Int_t sCharge = 0;
01218 if (gTrack->geometry()->charge() < 0) sCharge = 1;
01219
01220
01221
01222 StDedxPidTraits *pid, *pid70 = 0, *pidF = 0, *pid70U = 0, *pidFU = 0, *pid70A = 0;
01223 StProbPidTraits *pidprob = 0, *p = 0;
01224 Int_t NoFitPoints = gTrack->fitTraits().numberOfFitPoints(kTpcId);
01225 #ifdef CompareWithToF
01226 StBTofPidTraits* pidTof = 0;
01227 static const Int_t IdxH[4] = {kPidProton,kPidKaon,kPidPion,kPidElectron};
01228 #endif
01229 Double_t I70 = 0, D70 = 0, I70U = 0, D70U = 0, I70A = 0, D70A = 0;
01230 Double_t fitZ = 0, fitdZ = 1e10, fitZU = 0, fitdZU = 1e10;
01231 Int_t N70 = 0, NF = 0, N70A = 0;
01232 StSPtrVecTrackPidTraits &traits = gTrack->pidTraits();
01233 Double_t TrackLength70 = 0, TrackLength = 0;
01234 Double_t log2dX70 = 0, log2dX = 0;
01235 UInt_t size = traits.size();
01236 if (size) {
01237 for (UInt_t i = 0; i < size; i++) {
01238 if (! traits[i]) continue;
01239 if ( traits[i]->IsZombie()) continue;
01240 pid = dynamic_cast<StDedxPidTraits*>(traits[i]);
01241 if (pid && pid->detector() == kTpcId) {
01242 switch (pid->method()) {
01243 case kTruncatedMeanId:
01244 pid70 = pid; I70 = pid70->mean(); N70 = pid70->numberOfPoints();
01245 TrackLength70 = pid70->length(); D70 = pid70->errorOnMean();
01246 log2dX70 = pid70->log2dX();
01247 break;
01248 case kEnsembleTruncatedMeanId:
01249 pid70U = pid; I70U = pid70U->mean(); D70U = pid70U->errorOnMean();
01250 log2dX70 = pid70U->log2dX();
01251 break;
01252 case kLikelihoodFitId:
01253 pidF = pid;
01254 fitZ = TMath::Log(pidF->mean()); NF = pidF->numberOfPoints();
01255 TrackLength = pidF->length(); fitdZ = pidF->errorOnMean();
01256 log2dX = pidF->log2dX();
01257 break;
01258 case kWeightedTruncatedMeanId:
01259 pidFU = pid; fitZU = TMath::Log(pidFU->mean()); fitdZU = pidFU->errorOnMean();
01260 break;
01261 case kOtherMethodIdentifier:
01262 pid70A = pid; I70A = pid70A->mean(); D70A = pid70A->errorOnMean(); N70A = pid70A->numberOfPoints();
01263 break;
01264 default:
01265 break;
01266 }
01267 } else {
01268 p = dynamic_cast<StProbPidTraits*>(traits[i]);
01269 if (p) {pidprob = p; continue;}
01270 }
01271 }
01272 }
01273 #ifdef CompareWithToF
01274
01275 StTrackNode *node = gTrack->node();
01276 if (node) {
01277 StPrimaryTrack *pTrack = static_cast<StPrimaryTrack*>(node->track(primary));
01278 if (pTrack) {
01279 StSPtrVecTrackPidTraits &traitsp = pTrack->pidTraits();
01280 UInt_t size = traitsp.size();
01281 for (UInt_t i = 0; i < size; i++) {
01282 Short_t id = traitsp[i]->detector();
01283 if (id != kTofId) continue;
01284 StBTofPidTraits* p = dynamic_cast<StBTofPidTraits*>(traitsp[i]);
01285 if (p) {pidTof = p;}
01286 }
01287 }
01288 }
01289 #endif
01290 if (pid70 && ! pidF) {
01291 TrackLength = TrackLength70;
01292 log2dX = log2dX70;
01293 }
01294 Double_t PredB[NHYPS], Pred70B[NHYPS];
01295 Double_t PredBMN[2], Pred70BMN[2];
01296 Double_t date = GetDateTime().Convert() - timeOffSet;
01297 Double_t devZ[NHYPS], devZs[NHYPS], devToF[NHYPS];
01298 memset (devZ, 0, NHYPS*sizeof(Double_t));
01299 memset (devZs, 0, NHYPS*sizeof(Double_t));
01300 memset (devToF, 0, NHYPS*sizeof(Double_t));
01301 Double_t bg = TMath::Log10(pMomentum/StProbPidTraits::mPidParticleDefinitions[kPidPion]->mass());
01302 Double_t bghyp[NHYPS];
01303 Int_t l;
01304 PredBMN[0] = Pred70BMN[0] = 1;
01305 PredBMN[1] = Pred70BMN[1] = -1;
01306 for (l = kPidElectron; l < NHYPS; l += 1) {
01307 bghyp[l] = TMath::Log10(pMomentum*TMath::Abs(StProbPidTraits::mPidParticleDefinitions[l]->charge())/StProbPidTraits::mPidParticleDefinitions[l]->mass());
01308 PredB[l] = 1.e-6*StProbPidTraits::mPidParticleDefinitions[l]->charge()*StProbPidTraits::mPidParticleDefinitions[l]->charge()*
01309 TMath::Exp(m_Bichsel->GetMostProbableZ(bghyp[l]));
01310
01311 if (PredB[l] < PredBMN[0]) PredBMN[0] = PredB[l];
01312 if (PredB[l] > PredBMN[1]) PredBMN[1] = PredB[l];
01313 Pred70B[l] = 1.e-6*StProbPidTraits::mPidParticleDefinitions[l]->charge()*StProbPidTraits::mPidParticleDefinitions[l]->charge()*
01314 m_Bichsel->GetI70(bghyp[l]);
01315
01316 if (Pred70B[l] < Pred70BMN[0]) Pred70BMN[0] = Pred70B[l];
01317 if (Pred70B[l] > Pred70BMN[1]) Pred70BMN[1] = Pred70B[l];
01318
01319 if (pid70 && TrackLength70 > 40.) {
01320 hist70B[l][sCharge]->Fill(bghyp[l],TMath::Log(I70/Pred70B[l]));
01321 hitsB[l][sCharge]->Fill(bghyp[l],TMath::Log(PredB[l]));
01322 devZ[l] = TMath::Log(I70/Pred70B[l]);
01323 }
01324 if (pidF && TrackLength > 40.) {
01325 histzB[l][sCharge]->Fill(bghyp[l],fitZ - TMath::Log(PredB[l]));
01326 devZs[l] = TMath::Abs(devZ[l])/fitdZ;
01327 }
01328 }
01329 Double_t devZmin = 999;
01330 Int_t PiDStatus = 0;
01331 Int_t PiDkey = -1;
01332 Int_t PiDkeyU = -1;
01333 Int_t PiDkeyU3 = -1;
01334 Int_t lBest = -1;
01335 if (pidF && TrackLength > 40.) {
01336
01337 Double_t L10Mult = -1;
01338 if (St_trigDetSumsC::instance()) L10Mult = St_trigDetSumsC::instance()->mult();
01339
01340 StThreeVectorD xyz = gTrack->geometry()->helix().at(0);
01341 Double_t ZG = xyz.z();
01342 Double_t PhiDG = 180*xyz.phi();
01343 if (Pred70BMN[1] > 0 && D70 > 0) {
01344 if (TMath::Log(I70/Pred70BMN[0]) < -5*D70) {
01345 BaddEdxZPhi70[0]->Fill(ZG,PhiDG);
01346 BaddEdxMult70[0]->Fill(L10Mult);
01347 }
01348 if (TMath::Log(I70/Pred70BMN[1]) > 5*D70) {
01349 BaddEdxZPhi70[1]->Fill(ZG,PhiDG);
01350 BaddEdxMult70[1]->Fill(L10Mult);
01351 }
01352 }
01353 if (PredBMN[1] > 0 && fitdZ > 0) {
01354 if (fitZ - TMath::Log(PredBMN[0]) < -5*fitdZ) {
01355 BaddEdxZPhiZ[0]->Fill(ZG,PhiDG);
01356 BaddEdxMultZ[0]->Fill(L10Mult);
01357 }
01358 if (fitZ - TMath::Log(PredBMN[1]) > 5*fitdZ) {
01359 BaddEdxZPhiZ[1]->Fill(ZG,PhiDG);
01360 BaddEdxMultZ[1]->Fill(L10Mult);
01361 }
01362 }
01363 #ifdef CompareWithToF
01364
01365 if (pidTof) {
01366 for (l = kPidElectron; l < NHYPS; l++) {
01367 switch (l) {
01368 case kPidElectron:
01369 devToF[l] = pidTof->sigmaElectron();
01370 break;
01371 case kPidProton:
01372 devToF[l] = pidTof->sigmaProton();
01373 break;
01374 case kPidKaon:
01375 devToF[l] = pidTof->sigmaKaon();
01376 break;
01377 case kPidPion:
01378 devToF[l] = pidTof->sigmaPion();
01379 break;
01380 default:
01381 devToF[l] = 999.;
01382 break;
01383 }
01384 devToF[l] = TMath::Abs(devToF[l]);
01385 }
01386 for (l = kPidElectron; l < NHYPS; l++) {
01387 if (l == kPidMuon) continue;
01388 if (devToF[l] < 3.0) PiDStatus += 1<<l;
01389 if (devToF[l] < devZmin) {devZmin = devToF[l]; lBest = l;}
01390 }
01391 if (lBest >=0) {
01392 if (devToF[lBest] < 3.0) {
01393 PiDkey = lBest;
01394 Int_t lNext = -1;
01395 devZmin = 999;
01396 for (l = kPidElectron; l < NHYPS; l += 1) {
01397 if (l == kPidMuon) continue;
01398 if (l == lBest) continue;
01399 if (devToF[l] < devZmin) {devZmin = devToF[l]; lNext = l;}
01400 }
01401 if (lNext >= 0) {
01402 if (devToF[lNext] > 3.) {
01403 PiDkeyU = PiDkey;
01404 if (devToF[lNext] > 5.) PiDkeyU3 = PiDkeyU;
01405
01406 }
01407 }
01408 }
01409 }
01410 }
01411 #else
01412 for (l = kPidElectron; l < NHYPS; l++) {
01413 if (l == kPidMuon) continue;
01414 if (devZs[l] < 3.0) PiDStatus += 1<<l;
01415 if (devZs[l] < devZmin) {devZmin = devZs[l]; lBest = l;}
01416 }
01417 if (lBest >=0) {
01418 if (devZs[lBest] < 3.0) {
01419 PiDkey = lBest;
01420 Int_t lNext = -1;
01421 devZmin = 999;
01422 for (l = kPidElectron; l < NHYPS; l += 1) {
01423 if (l == kPidMuon) continue;
01424 if (l == lBest) continue;
01425 if (devZs[l] < devZmin) {devZmin = devZs[l]; lNext = l;}
01426 }
01427 if (lNext >= 0) {
01428 if (devZs[lNext] > 3.) {
01429 PiDkeyU = PiDkey;
01430 if (devZs[lNext] > 5.) PiDkeyU3 = PiDkeyU;
01431
01432 }
01433 }
01434 }
01435 }
01436 #endif
01437 }
01438 if (PiDkeyU3 >= 0) {
01439 l = PiDkeyU3;
01440 if (pid70 && TrackLength70 > 40.) {
01441 hist70BT[l][sCharge]->Fill(bghyp[l],TMath::Log(I70/Pred70B[l]));
01442 }
01443 if (pidF && TrackLength > 40.) {
01444 histzBT[l][sCharge]->Fill(bghyp[l],fitZ - TMath::Log(PredB[l]));
01445 }
01446 }
01447 if ((TESTBIT(m_Mode, kProbabilityPlot))) {
01448 if (TrackLength70 > 40) {
01449 Double_t Z70 = TMath::Log(I70/Pred70B[kPidPion]);
01450 Prob->Fill(bghyp[kPidPion],-0.05,Z70);
01451 if (pidprob) {
01452 Int_t N = pidprob->GetPidArray()->GetSize();
01453 for (int i = 0; i < N; i++) {
01454 Double_t p = pidprob->GetProbability(i);
01455 if (p > 0.99) p = 0.99;
01456 Prob->Fill(bghyp[kPidPion],i+p,Z70);
01457 }
01458 }
01459 }
01460 }
01461
01462 if (pidF && fitdZ > 0.05 && fitdZ < 0.15) {
01463 ffitZA->Fill(bg,devZ[kPidPion]);
01464 for (l = kPidElectron; l < NHYPS; l += 1) {
01465 ffitZ[l]->Fill(bg,devZ[kPidPion]);
01466 if (pidprob && pidprob->GetSum() > 0 && pidprob->GetProbability(l) > 0.9) ffitP[l]->Fill(bg,devZ[kPidPion]);
01467 if (PiDkeyU >= 0) ffitZU->Fill(bghyp[kPidPion],devZ[kPidPion]);
01468 if (PiDkeyU3 >= 0) ffitZU3->Fill(bghyp[kPidPion],devZ[kPidPion]);
01469 }
01470 }
01471 if (PredB[kPidPion] <= 0) {
01472 gMessMgr->Warning() << "StdEdxY2Maker:: Prediction for p = "
01473 << pMomentum << " and TrackLength = " << TrackLength
01474 << " is wrong = " << PredB[kPidPion] << " <<<<<<<<<<<<<" << endl;
01475 return;
01476 };
01477 if (pidF) {
01478 TPoints[0][0]->Fill(TrackLength,fitZ-TMath::Log(PredB[kPidPion]));
01479 #if 0
01480 Points[0][0]->Fill(NdEdx,(fitZ-TMath::Log(PredB[kPidPion]))/fitdZ);
01481 if (pMomentum > pMomin && pMomentum < pMomax) MPoints[0][0]->Fill(TrackLength,fitZ-TMath::Log(PredB[kPidPion]));
01482 #endif
01483 FitPull->Fill(TrackLength,(fitZ - TMath::Log(PredB[kPidPion]))/fitdZ);
01484 if (pidFU) {
01485 TPoints[0][2]->Fill(TrackLength,fitZU-TMath::Log(PredB[kPidPion]));
01486 #if 0
01487 Points[0][2]->Fill(NdEdx,(fitZU-TMath::Log(PredB[kPidPion]))/fitdZ);
01488 if (pMomentum > pMomin && pMomentum < pMomax) MPoints[0][2]->Fill(TrackLength,fitZU-TMath::Log(PredB[kPidPion]));
01489 #endif
01490 }
01491 }
01492 if (pid70) {
01493 TPoints[0][1]->Fill(TrackLength,TMath::Log(I70/Pred70B[kPidPion]));
01494 #if 0
01495 Points[0][1]->Fill(N70,TMath::Log(I70/PredB[kPidPion])/D70);
01496 if (pMomentum > pMomin && pMomentum < pMomax) MPoints[0][1]->Fill(TrackLength,TMath::Log(I70/Pred70B[kPidPion]));
01497 #endif
01498 Pull70->Fill(TrackLength,TMath::Log(I70/Pred70B[kPidPion])/D70);
01499 if (pid70U) {
01500 TPoints[0][3]->Fill(TrackLength,TMath::Log(I70U/Pred70B[kPidPion]));
01501 #if 0
01502 Points[0][3]->Fill(N70,TMath::Log(I70U/PredB[kPidPion])/D70);
01503 if (pMomentum > pMomin && pMomentum < pMomax) MPoints[0][3]->Fill(TrackLength,TMath::Log(I70U/Pred70B[kPidPion]));
01504 #endif
01505 }
01506 if (pid70A) {
01507 TPoints[0][4]->Fill(TrackLength,TMath::Log(I70A/Pred70B[kPidPion]));
01508 #if 0
01509 Points[0][4]->Fill(N70A,TMath::Log(I70A/PredB[kPidPion])/D70);
01510 if (pMomentum > pMomin && pMomentum < pMomax) MPoints[0][4]->Fill(TrackLength,TMath::Log(I70A/Pred70B[kPidPion]));
01511 #endif
01512 }
01513 }
01514 if (TrackLength > 20) {
01515
01516 Int_t k;
01517 Double_t bgL10 = TMath::Log10(pMomentum/StProbPidTraits::mPidParticleDefinitions[kPidPion]->mass());
01518 for (k = 0; k < NdEdx; k++) {
01519 FdEdx[k].zP = m_Bichsel->GetMostProbableZ(bgL10,1.);
01520
01521 FdEdx[k].sigmaP = m_Bichsel->GetRmsZ(bgL10,1.);
01522
01523 Double_t predB = 1.e-6*TMath::Exp(FdEdx[k].zP);
01524 FdEdx[k].dEdxN = TMath::Log(FdEdx[k].dEdx /predB);
01525 for (Int_t l = 0; l <= StTpcdEdxCorrection::kTpcLast; l++) {
01526 if (FdEdx[k].C[l].dEdx > 0)
01527 FdEdx[k].C[l].dEdxN = TMath::Log(FdEdx[k].C[l].dEdx/predB);
01528 }
01529 if (FdEdx[k].row < 14) {
01530 hdEI->Fill(TMath::Log10(FdEdx[k].dE));
01531 hdEUI->Fill(TMath::Log10(FdEdx[k].C[StTpcdEdxCorrection::kUncorrected].dE));
01532 hdERI->Fill(TMath::Log10(FdEdx[k].C[StTpcdEdxCorrection::kAdcCorrection].dE));
01533 hdEPI->Fill(TMath::Log10(FdEdx[k].C[StTpcdEdxCorrection::ktpcPressure].dE));
01534 hdETI->Fill(TMath::Log10(FdEdx[k].C[StTpcdEdxCorrection::ktpcTime].dE));
01535 hdESI->Fill(TMath::Log10(FdEdx[k].C[StTpcdEdxCorrection::kTpcSecRowB].dE));
01536 hdEZI->Fill(TMath::Log10(FdEdx[k].C[StTpcdEdxCorrection::kzCorrection].dE));
01537 hdEMI->Fill(TMath::Log10(FdEdx[k].C[StTpcdEdxCorrection::kMultiplicity].dE));
01538 }
01539 else {
01540 hdEO->Fill(TMath::Log10(FdEdx[k].dE));
01541 hdEUO->Fill(TMath::Log10(FdEdx[k].C[StTpcdEdxCorrection::kUncorrected].dE));
01542 hdERO->Fill(TMath::Log10(FdEdx[k].C[StTpcdEdxCorrection::kAdcCorrection].dE));
01543 hdEPO->Fill(TMath::Log10(FdEdx[k].C[StTpcdEdxCorrection::ktpcPressure].dE));
01544 hdETO->Fill(TMath::Log10(FdEdx[k].C[StTpcdEdxCorrection::ktpcTime].dE));
01545 hdESO->Fill(TMath::Log10(FdEdx[k].C[StTpcdEdxCorrection::kTpcSecRowB].dE));
01546 hdEZO->Fill(TMath::Log10(FdEdx[k].C[StTpcdEdxCorrection::kzCorrection].dE));
01547 hdEMO->Fill(TMath::Log10(FdEdx[k].C[StTpcdEdxCorrection::kMultiplicity].dE));
01548 }
01549 if ((TESTBIT(m_Mode, kAdcHistos))) {
01550 if (PiDkeyU3 >= 0) {
01551 Double_t betaXgamma = pMomentum*TMath::Abs(StProbPidTraits::mPidParticleDefinitions[PiDkeyU3]->charge())/StProbPidTraits::mPidParticleDefinitions[PiDkeyU3]->mass();
01552 Double_t zA = m_Bichsel->GetMostProbableZ(TMath::Log10(betaXgamma),1.) + 2*TMath::Log(TMath::Abs(StProbPidTraits::mPidParticleDefinitions[PiDkeyU3]->charge()));
01553 Double_t PredA = 1.e-6*TMath::Exp(zA);
01554 Double_t PredE = PredA*FdEdx[k].dx;
01555 Double_t PredEL = TMath::Log10(PredE);
01556 if (FdEdx[k].row < 14) {
01557 if (AdcI) AdcI->Fill(PredEL,TMath::Log10(FdEdx[k].C[StTpcdEdxCorrection::kUncorrected].dE));
01558 if (AdcIC) AdcIC->Fill(PredEL,TMath::Log10(FdEdx[k].dE));
01559 }
01560 else {
01561 if (AdcO) AdcO->Fill(PredEL,TMath::Log10(FdEdx[k].C[StTpcdEdxCorrection::kUncorrected].dE));
01562 if (AdcOC) AdcOC->Fill(PredEL,TMath::Log10(FdEdx[k].dE));
01563 }
01564 Double_t PredEU = PredE;
01565 Double_t PredEUL = TMath::Log10(PredEU);
01566 Double_t PredEULN = TMath::Log(PredEU);
01567 Double_t PredEUkeV = 1.e6*PredE;
01568 if (tpcGas) {
01569 if (FdEdx[k].row < 14) {
01570 if (Adc3I) Adc3I->Fill(PredEUL,TMath::Log10(FdEdx[k].C[StTpcdEdxCorrection::kUncorrected].dE));
01571 if (Adc3IC) Adc3IC->Fill(PredEUL,TMath::Log10(FdEdx[k].dE));
01572 if (Adc3Ip && Adc3Ip[PiDkeyU3])
01573 Adc3Ip[PiDkeyU3]->Fill(PredEUL,TMath::Log10(FdEdx[k].C[StTpcdEdxCorrection::kUncorrected].dE));
01574 if (sCharge == 0)
01575 AdcIZP->Fill(PredEUkeV,FdEdx[k].ZdriftDistanceO2,TMath::Log(FdEdx[k].C[StTpcdEdxCorrection::kUncorrected].dE)-PredEULN);
01576 else
01577 AdcIZN->Fill(PredEUkeV,FdEdx[k].ZdriftDistanceO2,TMath::Log(FdEdx[k].C[StTpcdEdxCorrection::kUncorrected].dE)-PredEULN);
01578 }
01579 else {
01580 if (Adc3O) Adc3O->Fill(PredEUL,TMath::Log10(FdEdx[k].C[StTpcdEdxCorrection::kUncorrected].dE));
01581 if (Adc3OC) Adc3OC->Fill(PredEUL,TMath::Log10(FdEdx[k].dE));
01582 if (Adc3Op && Adc3Op[PiDkeyU3])
01583 Adc3Op[PiDkeyU3]->Fill(PredEUL,TMath::Log10(FdEdx[k].C[StTpcdEdxCorrection::kUncorrected].dE));
01584 if (sCharge == 0)
01585 AdcOZP->Fill(PredEUkeV,FdEdx[k].ZdriftDistanceO2,TMath::Log(FdEdx[k].C[StTpcdEdxCorrection::kUncorrected].dE)-PredEULN);
01586 else
01587 AdcOZN->Fill(PredEUkeV,FdEdx[k].ZdriftDistanceO2,TMath::Log(FdEdx[k].C[StTpcdEdxCorrection::kUncorrected].dE)-PredEULN);
01588 }
01589 }
01590 }
01591 }
01592 if (pMomentum > pMomin && pMomentum < pMomax &&TrackLength > 40 ) {
01593 if (tpcGas) {
01594 if (inputTPCGasPressureP) inputTPCGasPressureP->Fill(tpcGas->inputTPCGasPressure,FdEdx[k].dEdxN);
01595 if (nitrogenPressureP) nitrogenPressureP->Fill(tpcGas->nitrogenPressure,FdEdx[k].dEdxN);
01596 if (gasPressureDiffP) gasPressureDiffP->Fill(tpcGas->gasPressureDiff,FdEdx[k].dEdxN);
01597 if (inputGasTemperatureP) inputGasTemperatureP->Fill(tpcGas->inputGasTemperature,FdEdx[k].dEdxN);
01598 if (outputGasTemperatureP) outputGasTemperatureP->Fill(tpcGas->outputGasTemperature,FdEdx[k].dEdxN);
01599 if (flowRateArgon1P) flowRateArgon1P->Fill(tpcGas->flowRateArgon1,FdEdx[k].dEdxN);
01600 if (flowRateArgon2P) flowRateArgon2P->Fill(tpcGas->flowRateArgon2,FdEdx[k].dEdxN);
01601 if (flowRateMethaneP) flowRateMethaneP->Fill(tpcGas->flowRateMethane,FdEdx[k].dEdxN);
01602 if (percentMethaneInP)
01603 percentMethaneInP->Fill(tpcGas->percentMethaneIn*1000./tpcGas->barometricPressure,
01604 FdEdx[k].C[StTpcdEdxCorrection::ktpcPressure-1].dEdxN);
01605 if (percentMethaneInPA)
01606 percentMethaneInPA->Fill(tpcGas->percentMethaneIn*1000./tpcGas->barometricPressure,
01607 FdEdx[k].C[StTpcdEdxCorrection::ktpcPressure].dEdxN);
01608 if (percentMethaneInPC)
01609 percentMethaneInPC->Fill(tpcGas->percentMethaneIn*1000./tpcGas->barometricPressure,
01610 FdEdx[k].dEdxN);
01611 if (ppmOxygenInP) ppmOxygenInP->Fill(tpcGas->ppmOxygenIn,FdEdx[k].dEdxN);
01612 if (flowRateExhaustP) flowRateExhaustP->Fill(tpcGas->flowRateExhaust,FdEdx[k].dEdxN);
01613 if (ppmWaterOutP) ppmWaterOutP->Fill(tpcGas->ppmWaterOut,FdEdx[k].C[StTpcdEdxCorrection::ktpcPressure-1].dEdxN);
01614 if (ppmWaterOutPA) ppmWaterOutPA->Fill(tpcGas->ppmWaterOut,FdEdx[k].C[StTpcdEdxCorrection::ktpcPressure].dEdxN);
01615 if (ppmWaterOutPC) ppmWaterOutPC->Fill(tpcGas->ppmWaterOut,FdEdx[k].dEdxN);
01616
01617 if (flowRateRecirculationP) flowRateRecirculationP->Fill(tpcGas->flowRateRecirculation,FdEdx[k].dEdxN);
01618 }
01619 if (St_trigDetSumsC::instance()) {
01620 if (St_trigDetSumsC::instance()->zdcX() > 0 && ZdcCP) ZdcCP->Fill(TMath::Log10(St_trigDetSumsC::instance()->zdcX()), FdEdx[k].dEdxN);
01621 if (St_trigDetSumsC::instance()->bbcYellowBkg() > 0 && bbcYellowBkg)
01622 bbcYellowBkg->Fill(TMath::Log10(St_trigDetSumsC::instance()->bbcYellowBkg()), FdEdx[k].dEdxN);
01623 if (St_trigDetSumsC::instance()->bbcBlueBkg() > 0 && bbcBlueBkg)
01624 bbcBlueBkg->Fill(TMath::Log10(St_trigDetSumsC::instance()->bbcBlueBkg()), FdEdx[k].dEdxN);
01625 if (St_trigDetSumsC::instance()->bbcX() > 0) {
01626 if (BBCP) BBCP->Fill(TMath::Log10(St_trigDetSumsC::instance()->bbcX()), FdEdx[k].dEdxN);
01627 if (BBC3) BBC3->Fill(FdEdx[k].row,TMath::Log10(St_trigDetSumsC::instance()->bbcX()), FdEdx[k].dEdxN);
01628 }
01629 if (St_trigDetSumsC::instance()->mult() > 0) {
01630 if (MultiplicityPI && FdEdx[k].row < 14) MultiplicityPI->Fill(TMath::Log10(St_trigDetSumsC::instance()->mult()), FdEdx[k].dEdxN);
01631 if (MultiplicityPO && FdEdx[k].row > 13) MultiplicityPO->Fill(TMath::Log10(St_trigDetSumsC::instance()->mult()), FdEdx[k].dEdxN);
01632 }
01633 }
01634 Double_t Pad2Edge = FdEdx[k].edge;
01635 if (TMath::Abs(Pad2Edge) > 5) {
01636 if (SecRow3 ) SecRow3->Fill(FdEdx[k].sector,FdEdx[k].row,FdEdx[k].C[StTpcdEdxCorrection::kTpcSecRowB-1].dEdxN);
01637 #if 0
01638 if (SecRow3A) SecRow3A->Fill(FdEdx[k].sector,FdEdx[k].row,FdEdx[k].C[StTpcdEdxCorrection::kTpcSecRowB].dEdxN);
01639 #endif
01640 if (SecRow3C) SecRow3C->Fill(FdEdx[k].sector,FdEdx[k].row,FdEdx[k].dEdxN);
01641 }
01642 if (Zdc3C && FdEdx[k].Zdc > 0) Zdc3C->Fill(FdEdx[k].row,TMath::Log10(FdEdx[k].Zdc),FdEdx[k].dEdxN);
01643
01644 Double_t xyzD[3] = {FdEdx[k].xyzD[0],FdEdx[k].xyzD[1],FdEdx[k].xyzD[2]};
01645
01646 Double_t PhiD = 180./TMath::Pi()*TMath::ATan2(xyzD[0],xyzD[1]);
01647 Double_t ThetaD = 180./TMath::Pi()*TMath::ATan2(-xyzD[2],TMath::Sqrt(xyzD[0]*xyzD[0]+xyzD[1]*xyzD[1]));
01648 if (Phi3) Phi3->Fill(FdEdx[k].PhiR,FdEdx[k].row,FdEdx[k].dEdxN);
01649 if (Phi3D) Phi3D->Fill(PhiD,FdEdx[k].row,FdEdx[k].dEdxN);
01650 if (Theta3D) Theta3D->Fill(ThetaD,FdEdx[k].row,FdEdx[k].dEdxN);
01651 if ((TESTBIT(m_Mode, kMip))) {
01652 if (SecRow3Mip && TMath::Abs(devZs[kPidPion]) < 2)
01653 SecRow3Mip->Fill(FdEdx[k].row,
01654 TMath::Log(FdEdx[k].dx)/TMath::Log(2.),
01655 FdEdx[k].dEdxN);
01656 }
01657 if (tpcGas && Pressure) {
01658 Double_t p = tpcGas->barometricPressure;
01659 Double_t t = tpcGas->inputGasTemperature/298.2;
01660 if (p > 0) {
01661 Double_t press = TMath::Log(p);
01662 if (Pressure) Pressure ->Fill(FdEdx[k].row,press,FdEdx[k].C[StTpcdEdxCorrection::ktpcPressure-1].dEdxN);
01663 if (PressureA) PressureA->Fill(FdEdx[k].row,press,FdEdx[k].C[StTpcdEdxCorrection::ktpcPressure].dEdxN);
01664 if (PressureC) PressureC->Fill(FdEdx[k].row,press,FdEdx[k].dEdxN);
01665 }
01666 Double_t V = St_tpcAnodeHVavgC::instance()->voltagePadrow(FdEdx[k].sector,FdEdx[k].row);
01667 Double_t VN = (FdEdx[k].row <= 13) ? V - 1170 : V - 1390;
01668 if (V > 0) {
01669 if (Voltage) Voltage ->Fill(45*(FdEdx[k].sector-1)+FdEdx[k].row,
01670 VN,FdEdx[k].C[StTpcdEdxCorrection::ktpcPressure-1].dEdxN);
01671 if (VoltageC) VoltageC->Fill(FdEdx[k].row,
01672 V,FdEdx[k].dEdxN);
01673 }
01674 Double_t Q = St_tpcAvCurrentC::instance()->chargeI();
01675 if (FdEdx[k].row > 13) Q = St_tpcAvCurrentC::instance()->chargeO();
01676 AcCharge->Fill(FdEdx[k].row,Q,FdEdx[k].dEdxN);
01677 Qcm->Fill(FdEdx[k].row,FdEdx[k].Qcm,FdEdx[k].dEdxN);
01678 AvCurrent->Fill(FdEdx[k].row,FdEdx[k].Crow,FdEdx[k].dEdxN);
01679 if (p*t > 0) {
01680 Double_t temp = TMath::Log(p/t);
01681 if (PressureT) PressureT ->Fill(FdEdx[k].row,temp,FdEdx[k].C[StTpcdEdxCorrection::ktpcPressure-1].dEdxN);
01682 if (PressureTA) PressureTA->Fill(FdEdx[k].row,temp,FdEdx[k].C[StTpcdEdxCorrection::ktpcPressure].dEdxN);
01683 if (PressureTC) PressureTC->Fill(FdEdx[k].row,temp,FdEdx[k].dEdxN);
01684 }
01685 }
01686 if (Time) Time->Fill(date,FdEdx[k].C[StTpcdEdxCorrection::ktpcTime-1].dEdxN);
01687 if (TimeP) TimeP->Fill(date,FdEdx[k].C[StTpcdEdxCorrection::ktpcTime].dEdxN);
01688 if (TimeC) TimeC->Fill(date,FdEdx[k].dEdxN);
01689 #if 0
01690 if (Z3O) Z3O->Fill(FdEdx[k].row,FdEdx[k].ZdriftDistanceO2,FdEdx[k].C[StTpcdEdxCorrection::kDrift-1].dEdxN);
01691 if (Z3OA)Z3OA->Fill(FdEdx[k].row,FdEdx[k].ZdriftDistanceO2,FdEdx[k].C[StTpcdEdxCorrection::kDrift].dEdxN);
01692 #endif
01693 if (Z3OC)Z3OC->Fill(FdEdx[k].row,FdEdx[k].ZdriftDistanceO2,FdEdx[k].dEdxN);
01694 #if 0
01695 if (Z3OW) Z3OW->Fill(FdEdx[k].row,FdEdx[k].ZdriftDistanceO2W,FdEdx[k].C[StTpcdEdxCorrection::kDrift-1].dEdxN);
01696 if (Z3OWA)Z3OWA->Fill(FdEdx[k].row,FdEdx[k].ZdriftDistanceO2W,FdEdx[k].C[StTpcdEdxCorrection::kDrift].dEdxN);
01697 if (Z3OWC)Z3OWC->Fill(FdEdx[k].row,FdEdx[k].ZdriftDistanceO2W,FdEdx[k].dEdxN);
01698 if (Z3) Z3->Fill(FdEdx[k].row,FdEdx[k].ZdriftDistance, FdEdx[k].C[StTpcdEdxCorrection::kzCorrection-1].dEdxN);
01699 if (Z3A) Z3A->Fill(FdEdx[k].row,FdEdx[k].ZdriftDistance, FdEdx[k].C[StTpcdEdxCorrection::kzCorrection].dEdxN);
01700 #endif
01701 if (Z3C) Z3C->Fill(FdEdx[k].row,FdEdx[k].ZdriftDistance, FdEdx[k].dEdxN);
01702 #if 0
01703 if (ETA) ETA->Fill(FdEdx[k].row,eta,FdEdx[k].dEdxN);
01704 if (ETA3)ETA3->Fill(FdEdx[k].row,eta,FdEdx[k].dEdxN);
01705 #endif
01706 if (Edge3) Edge3->Fill(FdEdx[k].row,FdEdx[k].edge, FdEdx[k].dEdxN);
01707 if (St_trigDetSumsC::instance() && St_trigDetSumsC::instance()->mult() > 0) {
01708 if (MulRow) MulRow->Fill(TMath::Log10(St_trigDetSumsC::instance()->mult()),FdEdx[k].row,FdEdx[k].C[StTpcdEdxCorrection::kMultiplicity-1].dEdxN);
01709 if (MulRowC) MulRowC->Fill(TMath::Log10(St_trigDetSumsC::instance()->mult()),FdEdx[k].row,FdEdx[k].dEdxN);
01710 }
01711 if (dXdE ) dXdE->Fill(FdEdx[k].dx,FdEdx[k].row,FdEdx[k].C[StTpcdEdxCorrection::kdXCorrection-1].dEdxN);
01712 if (dXdEA) dXdEA->Fill(FdEdx[k].dx,FdEdx[k].row,FdEdx[k].C[StTpcdEdxCorrection::kdXCorrection].dEdxN);
01713 if (dXdEC) dXdEC->Fill(FdEdx[k].dx,FdEdx[k].row,FdEdx[k].dEdxN);
01714 }
01715
01716 if (TESTBIT(m_Mode, kZBGX) && zbgx)
01717 zbgx[2*PiDkeyU3+sCharge]->Fill(bghyp[PiDkeyU3],TMath::Log(FdEdx[k].dx)/TMath::Log(2.),FdEdx[k].dEdxL-GeV2keV);
01718 }
01719 if ((TESTBIT(m_Mode, kCORRELATION))) Correlations();
01720 }
01721 #ifdef CompareWithToF
01722
01723 if (NdEdx > 0 && ftrack && ftree && TrackLength70 > 20.0 && pidTof) {
01724 ftrack->Clear();
01725 ftrack->sCharge = 1 - 2*sCharge;
01726 ftrack->p = pMomentum;
01727 ftrack->pX = g3.x();
01728 ftrack->pY = g3.y();
01729 ftrack->pZ = g3.z();
01730 ftrack->R0 = gTrack->geometry()->origin().mag();
01731 ftrack->Z0 = gTrack->geometry()->origin().z();
01732 ftrack->Phi0 = gTrack->geometry()->origin().phi();
01733 ftrack->NoFitPoints = NoFitPoints;
01734 ftrack->NdEdx = NdEdx;
01735 ftrack->N70 = N70;
01736 ftrack->I70 = I70;
01737 ftrack->D70 = D70;
01738 ftrack->TrackLength70 = TrackLength70;
01739 ftrack->fitZ = fitZ;
01740 ftrack->fitdZ = fitdZ;
01741 ftrack->TrackLength = TrackLength;
01742 Double_t *h = 0;
01743 for (Int_t l = 0; l < 4; l++) {
01744 Int_t k = IdxH[l];
01745 if (k == kPidProton ) h = &ftrack->PredP;
01746 else if (k == kPidKaon ) h = &ftrack->PredK;
01747 else if (k == kPidPion ) h = &ftrack->Predpi;
01748 else if (k == kPidElectron) h = &ftrack->PredE;
01749 else continue;
01750 h[0] = PredB[k];
01751 h[1] = Pred70B[k];
01752 h[2] = devZs[k];
01753 h[3] = devZ[k];
01754 h[4] = devToF[k];
01755 if (pidprob) h[5] = pidprob->GetProbability(k);
01756 else h[5] = -1;
01757 }
01758 for (Int_t k = 0; k < NdEdx; k++) {
01759 ftrack->AddPoint(FdEdx[k]);
01760 }
01761 ftree->Fill();
01762 }
01763 #endif
01764 return;
01765 }
01766
01767 void StdEdxY2Maker::PrintdEdx(Int_t iop) {
01768 const Int_t NOpts = 20;
01769 const Char_t *Names[NOpts] = {"CdEdx","FdEdx","dEdxS","dEdxU","dEdxR",
01770 "dEdxS","dEdxS","dEdxP","dEdxt","dEdxO",
01771 "dEdxM","dEdxZ","dEdxm","dEdxT","dEdxW",
01772 "dEdxC","dEdxE","dEdxp","dEdxX","dEdxd"};
01773 if (iop < 0 || iop >= NOpts) return;
01774 dEdxY2_t *pdEdx = 0;
01775 Double_t dEdx;
01776 Double_t I = 0;
01777 Int_t N70 = NdEdx - (int) (0.3*NdEdx + 0.5);
01778 Int_t N60 = NdEdx - (int) (0.4*NdEdx + 0.5);
01779 Double_t I70 = 0, I60 = 0;
01780 Double_t avrz = 0;
01781 for (int i=0; i< NdEdx; i++) {
01782 dEdx = 0;
01783 if (iop == 0) {pdEdx = &CdEdx[i]; dEdx = CdEdx[i].dEdx;}
01784 else if (iop == 1) {pdEdx = &FdEdx[i]; dEdx = FdEdx[i].dEdx;}
01785 else if (iop == 2) {pdEdx = &dEdxS[i]; dEdx = dEdxS[i].dEdx;}
01786 else if (iop >= 3) {pdEdx = &FdEdx[i]; dEdx = FdEdx[i].C[StTpcdEdxCorrection::kUncorrected+iop-3].dEdx;}
01787 I = (i*I + pdEdx->dEdx)/(i+1);
01788
01789
01790
01791
01792
01793
01794
01795
01796 cout << Form("%s %2i S/R %2i/%2i dEdx(keV/cm) %8.2f dx %5.2f x[%8.2f,%8.2f,%8.2f]",
01797 Names[iop],i,pdEdx->sector,pdEdx->row,1.e6*dEdx, pdEdx->dx, pdEdx->xyz[0], pdEdx->xyz[1], pdEdx->xyz[2]);
01798 cout << Form(" d[%8.2f,%8.2f,%8.2f] Sum %8.2f Prob %8.5f", pdEdx->xyzD[0], pdEdx->xyzD[1], pdEdx->xyzD[2],1.e6*I,pdEdx->Prob) << endl;
01799 if (iop == 2) {
01800 if (i < N60) I60 += dEdx;
01801 if (i < N70) I70 += dEdx;
01802 if (i == N60 - 1) {
01803 I60 /= N60;
01804 cout << " ======================= I60 \t" << I60 << endl;
01805 }
01806 if (i == N70 - 1) {
01807 I70 /= N70;
01808 cout << " ======================= I70 \t" << I70 << endl;
01809 }
01810 }
01811 avrz += TMath::Log(dEdx);
01812 }
01813 if (NdEdx) avrz /= NdEdx;
01814 cout << "mean dEdx \t" << I << "\tExp(avrz)\t" << TMath::Exp(avrz) << endl;
01815 }
01816
01817 Double_t StdEdxY2Maker::LikeliHood(Double_t Xlog10bg, Int_t NdEdx, dEdxY2_t *dEdx, Double_t chargeSq) {
01818
01819
01820 const static Double_t ProbCut = 1.e-4;
01821 const static Double_t GeV2keV = TMath::Log(1.e-6);
01822 Double_t f = 0;
01823 #if 0
01824 St_tpcCorrectionC *dXCorrection = 0;
01825 if (m_TpcdEdxCorrection) dXCorrection = m_TpcdEdxCorrection->dXCorrection();
01826 Int_t k = 0;
01827 if (dXCorrection) {
01828 tpcCorrection_st *dXC = ((St_tpcCorrection *)dXCorrection->Table())->GetTable();
01829 if (dXC) {
01830 Int_t N = dXC->nrows;
01831 if (N > 3) k = 3;
01832 }
01833 }
01834 #endif
01835 for (int i=0;i<NdEdx; i++) {
01836
01837 Double_t Ylog2dx = 1;
01838 Double_t sigmaC = 0;
01839 #if 0
01840 StTpcdEdxCorrection::ESector l = StTpcdEdxCorrection::kTpcInner;
01841 if (dEdx[i].row > NumberOfInnerRows) l = StTpcdEdxCorrection::kTpcOuter;
01842 if (dXCorrection) sigmaC = dXCorrection->CalcCorrection(k,Ylog2dx);
01843 #endif
01844 Double_t zMostProb = m_Bichsel->GetMostProbableZ(Xlog10bg,Ylog2dx) + TMath::Log(chargeSq);
01845 Double_t sigma = m_Bichsel->GetRmsZ(Xlog10bg,Ylog2dx) + sigmaC;
01846 Double_t xi = (dEdx[i].dEdxL - GeV2keV - zMostProb)/sigma;
01847 Double_t Phi = m_Bichsel->GetProbability(Xlog10bg,Ylog2dx,xi);
01848 dEdx[i].Prob = Phi/sigma;
01849 if (dEdx[i].Prob < ProbCut) {
01850 dEdx[i].Prob = ProbCut;
01851 }
01852 f -= 2*TMath::Log( dEdx[i].Prob );
01853 }
01854 return f;
01855 }
01856
01857 void StdEdxY2Maker::Landau(Double_t x, Double_t *val){
01858
01859 Double_t params[9] = {
01860 -7.74975e+00,
01861 6.53414e+00,
01862 1.21524e+00,
01863 3.31409e+00,
01864 -2.58291e+00,
01865 3.51463e+00,
01866 -3.47755e+00,
01867 3.77698e-02,
01868 6.67913e-01};
01869 Double_t dev1 = (x-params[1])/params[2];
01870 Double_t dev2 = (x-params[4])/params[5];
01871 Double_t dev3 = (x-params[7])/params[8];
01872 Double_t d = TMath::Exp(params[6]-0.5*dev3*dev3);
01873 Double_t dp = -dev3/params[8]*d;
01874 Double_t c = TMath::Exp(params[3]-0.5*dev2*dev2 + d);
01875 Double_t cp = (-dev2/params[5] + dp)*c;
01876 val[0] = params[0]-0.5*dev1*dev1+c;
01877 val[1] = - dev1/params[2]+cp;
01878 }
01879
01880 void StdEdxY2Maker::fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag) {
01881 Double_t Val[2];
01882
01883
01884 static Double_t sigma_p[3] = { 5.66734e-01, -1.24725e-01, 1.96085e-02};
01885 f = 0.;
01886 gin[0] = 0.;
01887 gin[1] = 0.;
01888 #if 0
01889 Double_t weightT = 0;
01890 #endif
01891 for (int i=0;i<NdEdx; i++) {
01892
01893 Double_t X = TMath::Log(FdEdx[i].dx);
01894 Double_t sigma = sigma_p[2];
01895 for (int n = 1; n>=0; n--) sigma = X*sigma + sigma_p[n];
01896 FdEdx[i].zdev = (FdEdx[i].dEdxL-par[0])/sigma;
01897 Landau(FdEdx[i].zdev,Val);
01898 FdEdx[i].Prob = TMath::Exp(Val[0]);
01899 #if 0
01900 Double_t weight = FdEdx[i].Weight;
01901 weightT += weight;
01902 f -= Val[0]*weight;
01903 gin[0] += Val[1]/sigma*weight;
01904 #else
01905 f -= Val[0];
01906 gin[0] += Val[1]/sigma;
01907 #endif
01908 }
01909 #if 0
01910 weightT /= NdEdx;
01911 f /= weightT;
01912 gin[0] /= weightT;
01913 #endif
01914 }
01915
01916 void StdEdxY2Maker::DoFitZ(Double_t &chisq, Double_t &fitZ, Double_t &fitdZ){
01917 Double_t avz = 0;
01918 for (int i=0;i<NdEdx;i++) avz += FdEdx[i].dEdxL;
01919 if (NdEdx>5) {
01920 avz /= NdEdx;
01921 Double_t arglist[10];
01922 Int_t ierflg = 0;
01923 m_Minuit->SetFCN(fcn);
01924
01925 if (Debug() < 2) {
01926 arglist[0] = -1;
01927 m_Minuit->mnexcm("set print",arglist, 1, ierflg);
01928 }
01929 m_Minuit->mnexcm("set NOW",arglist, 0, ierflg);
01930 m_Minuit->mnexcm("CLEAR",arglist, 0, ierflg);
01931 arglist[0] = 0.5;
01932 m_Minuit->mnexcm("SET ERR", arglist ,1,ierflg);
01933 m_Minuit->mnparm(0, "mean", avz, 0.01,0.,0.,ierflg);
01934 if (Debug() < 2) arglist[0] = 1.;
01935 else arglist[0] = 0.;
01936 m_Minuit->mnexcm("SET GRAD",arglist,1,ierflg);
01937 arglist[0] = 500;
01938 arglist[1] = 1.;
01939 m_Minuit->mnexcm("MIGRAD", arglist ,2,ierflg);
01940 m_Minuit->mnexcm("HESSE ",arglist,0,ierflg);
01941 Double_t edm,errdef;
01942 Int_t nvpar,nparx,icstat;
01943 m_Minuit->mnstat(chisq,edm,errdef,nvpar,nparx,icstat);
01944 m_Minuit->GetParameter(0, fitZ, fitdZ);
01945 }
01946 else {
01947 fitZ = fitdZ = chisq = -999.;
01948 }
01949 }
01950
01951 void StdEdxY2Maker::TrigHistos(Int_t iok) {
01952 static TProfile *BarPressure = 0, *inputTPCGasPressure = 0;
01953 static TProfile *nitrogenPressure = 0, *gasPressureDiff = 0, *inputGasTemperature = 0;
01954 static TProfile *outputGasTemperature = 0, *flowRateArgon1 = 0, *flowRateArgon2 = 0, *flowRateMethane = 0;
01955 static TProfile *percentMethaneIn = 0, *ppmOxygenIn = 0, *flowRateExhaust = 0;
01956 static TProfile *ppmWaterOut = 0, *ppmOxygenOut = 0, *flowRateRecirculation = 0;
01957
01958 static TH1F *Multiplicity;
01959 #if 0
01960 static TH2F *Zdc = 0;
01961 #endif
01962 static TH1F *ZdcC = 0;
01963 static TH1F *BBC = 0;
01964 if (! iok && !BarPressure) {
01965 TDatime t1(tMin,0);
01966 TDatime t2(tMax,0);
01967 UInt_t i1 = t1.Convert() - timeOffSet;
01968 UInt_t i2 = t2.Convert() - timeOffSet;
01969 Int_t Nt = (i2 - i1)/(3600);
01970 BarPressure = new TProfile("BarPressure","barometricPressure (mbar) versus time",Nt,i1,i2);
01971 inputTPCGasPressure = new TProfile("inputTPCGasPressure","inputTPCGasPressure (mbar) versus time",Nt,i1,i2);
01972 nitrogenPressure = new TProfile("nitrogenPressure","nitrogenPressure (mbar) versus time",Nt,i1,i2);
01973 gasPressureDiff = new TProfile("gasPressureDiff","gasPressureDiff (mbar) versus time",Nt,i1,i2);
01974 inputGasTemperature = new TProfile("inputGasTemperature","inputGasTemperature (degrees C) versus time",Nt,i1,i2);
01975 outputGasTemperature = new TProfile("outputGasTemperature","outputGasTemperature (degrees C) versus time",Nt,i1,i2);
01976 flowRateArgon1 = new TProfile("flowRateArgon1","flowRateArgon1 (liters/min) versus time",Nt,i1,i2);
01977 flowRateArgon2 = new TProfile("flowRateArgon2","flowRateArgon2 (liters/min) versus time",Nt,i1,i2);
01978 flowRateMethane = new TProfile("flowRateMethane","flowRateMethane (liters/min) versus time",Nt,i1,i2);
01979 percentMethaneIn = new TProfile("percentMethaneIn","percentMethaneIn (percent) versus time",Nt,i1,i2);
01980 ppmOxygenIn = new TProfile("ppmOxygenIn","ppmOxygenIn (ppm) versus time",Nt,i1,i2);
01981 flowRateExhaust = new TProfile("flowRateExhaust","flowRateExhaust (liters/min) versus time",Nt,i1,i2);
01982 ppmWaterOut = new TProfile("ppmWaterOut","ppmWaterOut (ppm) versus time",Nt,i1,i2);
01983 ppmOxygenOut = new TProfile("ppmOxygenOut","ppmOxygenOut (ppm) versus time",Nt,i1,i2);
01984 flowRateRecirculation = new TProfile("flowRateRecirculation","flowRateRecirculation (liters/min) versus time",Nt,i1,i2);
01985
01986
01987
01988
01989
01990 #if 0
01991 Zdc = new TH2F("Zdc","ZdcEastRate versus ZdcWestRate (log10)",100,0,10,100,0,10);
01992 #endif
01993 ZdcC = new TH1F("ZdcC","ZdcCoincidenceRate (log10)",100,0,10);
01994 Multiplicity = new TH1F("Multiplicity","Multiplicity (log10)",100,0,10);
01995 BBC = new TH1F("BBC","BbcCoincidenceRate (log10)",100,0,10);
01996 }
01997 else {
01998 tpcGas_st *tpcgas = m_TpcdEdxCorrection && m_TpcdEdxCorrection->tpcGas() ? m_TpcdEdxCorrection->tpcGas()->GetTable():0;
01999 Double_t date = GetDateTime().Convert() - timeOffSet;
02000 if (tpcgas) {
02001 if (BarPressure) BarPressure->Fill(date,tpcgas->barometricPressure);
02002 if (inputTPCGasPressure) inputTPCGasPressure->Fill(date,tpcgas->inputTPCGasPressure);
02003 if (nitrogenPressure) nitrogenPressure->Fill(date,tpcgas->nitrogenPressure);
02004 if (gasPressureDiff) gasPressureDiff->Fill(date,tpcgas->gasPressureDiff);
02005 if (inputGasTemperature) inputGasTemperature->Fill(date,tpcgas->inputGasTemperature);
02006 if (outputGasTemperature) outputGasTemperature->Fill(date,tpcgas->outputGasTemperature);
02007 if (flowRateArgon1) flowRateArgon1->Fill(date,tpcgas->flowRateArgon1);
02008 if (flowRateArgon2) flowRateArgon2->Fill(date,tpcgas->flowRateArgon2);
02009 if (flowRateMethane) flowRateMethane->Fill(date,tpcgas->flowRateMethane);
02010 if (percentMethaneIn)
02011 percentMethaneIn->Fill(date,tpcgas->percentMethaneIn*1000./tpcgas->barometricPressure);
02012 if (ppmOxygenIn) ppmOxygenIn->Fill(date,tpcgas->ppmOxygenIn);
02013 if (flowRateExhaust) flowRateExhaust->Fill(date,tpcgas->flowRateExhaust);
02014 if (ppmWaterOut) ppmWaterOut->Fill(date,tpcgas->ppmWaterOut);
02015 if (ppmOxygenOut) ppmOxygenOut->Fill(date,tpcgas->ppmOxygenOut);
02016 if (flowRateRecirculation) flowRateRecirculation->Fill(date,tpcgas->flowRateRecirculation);
02017 }
02018 if (St_trigDetSumsC::instance()) {
02019 #if 0
02020 if (St_trigDetSumsC::instance()->zdcWest() > 0 &&
02021 St_trigDetSumsC::instance()->zdcEast() > 0 && Zdc) Zdc->Fill(TMath::Log10(St_trigDetSumsC::instance()->zdcWest()),
02022 TMath::Log10(St_trigDetSumsC::instance()->zdcEast()));
02023 #endif
02024 if (St_trigDetSumsC::instance()->zdcX() > 0 && ZdcC) ZdcC->Fill(TMath::Log10(St_trigDetSumsC::instance()->zdcX()));
02025 if (St_trigDetSumsC::instance()->bbcX() > 0 && BBC) BBC->Fill(TMath::Log10(St_trigDetSumsC::instance()->bbcX()));
02026 if (St_trigDetSumsC::instance()->mult() > 0 && Multiplicity) Multiplicity->Fill(TMath::Log10(St_trigDetSumsC::instance()->mult()));
02027 }
02028 }
02029 }
02030
02031 void StdEdxY2Maker::XyzCheck(StGlobalCoordinate *global, Int_t iokCheck) {
02032 static TH3F *XYZ = 0, *XYZbad = 0;
02033 if (! global && !XYZ) {
02034 if (Debug()) gMessMgr->Warning() << "StdEdxY2Maker::XyzCheck XYZ check Histograms" << endm;
02035 XYZ = new TH3F("XYZ","xyz for clusters",80,-200,200,80,-200,200,84,-210,210);
02036 XYZbad = new TH3F("XYZbad","xyz for clusters with mismatched sectors",
02037 80,-200,200,80,-200,200,84,-210,210);
02038 }
02039 else
02040 if (XYZ) XYZ->Fill( global->position().x(), global->position().y(), global->position().z());
02041 if (iokCheck && XYZbad) XYZbad->Fill( global->position().x(), global->position().y(), global->position().z());
02042 }
02043
02044 void StdEdxY2Maker::QAPlots(StGlobalTrack* gTrack) {
02045 static TH2F *fTdEdxP70 = 0, *fTdEdxP70pi = 0, *fTdEdxP70e = 0, *fTdEdxP70K = 0, *fTdEdxP70P = 0;
02046 static TH2F *fTdEdxPF = 0, *fTdEdxPFpi = 0, *fTdEdxPFe = 0, *fTdEdxPFK = 0, *fTdEdxPFP = 0;
02047 static StTpcDedxPidAlgorithm PidAlgorithm;
02048 static StElectron* Electron = StElectron::instance();
02049 static StPionPlus* Pion = StPionPlus::instance();
02050 static StKaonPlus* Kaon = StKaonPlus::instance();
02051 static StProton* Proton = StProton::instance();
02052 static const Double_t Log10E = TMath::Log10(TMath::Exp(1.));
02053 static int first=0;
02054 if (! gTrack) {
02055 TFile *f = 0;
02056 if (TESTBIT(m_Mode, kCalibration)) {
02057 f = GetTFile();
02058 if (f) f->cd();
02059 }
02060 if (!first) {
02061 fZOfBadHits = new TH1F*[fNZOfBadHits];
02062 static const Char_t *BadCaseses[fNZOfBadHits] =
02063 {"it is not used in track fit",
02064 "it is flagged ",
02065 "track length is inf ",
02066 "it does not pass check ",
02067 "dx is out interval [0.5,25]",
02068 "Sector/Row gain < 0",
02069 "drift distance < min || drift distance > max",
02070 "dE < 0 or dx < 0",
02071 "Edge effect",
02072 "Anode Voltage problem",
02073 "Total no.of rejected clusters"
02074 };
02075 for (Int_t i = 0; i < fNZOfBadHits; i++)
02076 fZOfBadHits[i] = new TH1F(Form("ZOfBadHits%i",i),
02077 Form("Z of rejected clusters because %s",BadCaseses[i]),
02078 100,-210,210);
02079 fZOfGoodHits = new TH1F("ZOfGoodHits","Z of accepted clusters",100,-210,210);
02080 fPhiOfBadHits = new TH1F("PhiOfBadHits","Phi of rejected clusters",100, -TMath::Pi(), TMath::Pi());
02081 fTracklengthInTpcTotal = new TH1F("TracklengthInTpcTotal","Total track in TPC",100,0,200);
02082 fTracklengthInTpc = new TH1F("TracklengthInTpc","Track length in TPC used for dE/dx",100,0,200);
02083 fTdEdxPF = new TH2F("TdEdxPF","log10(dE/dx(fit)(keV/cm)) versus log10(p(GeV/c)) for Tpc TrackLength > 40 cm",
02084 250,-1.,4., 500,0.,2.5);
02085 fTdEdxP70 = new TH2F("TdEdxP70","log10(dE/dx(I70)(keV/cm)) versus log10(p(GeV/c)) for Tpc TrackLength > 40 cm",
02086 250,-1.,4., 500,0.,2.5);
02087 fTdEdxP70pi = new TH2F("TdEdxP70pi","log10(dE/dx(I70)(keV/cm)) versus log10(p(GeV/c)) for Tpc TrackLength > 40 cm |nSigmaPion| < 1",
02088 250,-1.,4., 500,0.,2.5);
02089 fTdEdxP70pi->SetMarkerColor(2);
02090 fTdEdxP70e = new TH2F("TdEdxP70e","log10(dE/dx(I70)(keV/cm)) versus log10(p(GeV/c)) for Tpc TrackLength > 40 cm |nSigmaElectron| < 1",
02091 250,-1.,4., 500,0.,2.5);
02092 fTdEdxP70e->SetMarkerColor(3);
02093 fTdEdxP70K = new TH2F("TdEdxP70K","log10(dE/dx(I70)(keV/cm)) versus log10(p(GeV/c)) for Tpc TrackLength > 40 cm |nSigmaKaon| < 1",
02094 250,-1.,4., 500,0.,2.5);
02095 fTdEdxP70K->SetMarkerColor(4);
02096 fTdEdxP70P = new TH2F("TdEdxP70P","log10(dE/dx(I70)(keV/cm)) versus log10(p(GeV/c)) for Tpc TrackLength > 40 cm |nSigmaProton| < 1",
02097 250,-1.,4., 500,0.,2.5);
02098 fTdEdxP70P->SetMarkerColor(6);
02099
02100 fTdEdxPFpi = new TH2F("TdEdxPFpi","log10(dE/dx(Ifit)(keV/cm)) versus log10(p(GeV/c)) for Tpc TrackLength > 40 cm and Prob > 68.3",
02101 250,-1.,4., 500,0.,2.5);
02102 fTdEdxPFpi->SetMarkerColor(2);
02103 fTdEdxPFe = new TH2F("TdEdxPFe","log10(dE/dx(Ifit)(keV/cm)) versus log10(p(GeV/c)) for Tpc TrackLength > 40 cm and Prob > 68.3",
02104 250,-1.,4., 500,0.,2.5);
02105 fTdEdxPFe->SetMarkerColor(3);
02106 fTdEdxPFK = new TH2F("TdEdxPFK","log10(dE/dx(Ifit)(keV/cm)) versus log10(p(GeV/c)) for Tpc TrackLength > 40 cm and Prob > 68.3",
02107 250,-1.,4., 500,0.,2.5);
02108 fTdEdxPFK->SetMarkerColor(4);
02109 fTdEdxPFP = new TH2F("TdEdxPFP","log10(dE/dx(Ifit)(keV/cm)) versus log10(p(GeV/c)) for Tpc TrackLength > 40 cm and Prob > 68.3",
02110 250,-1.,4., 500,0.,2.5);
02111 fTdEdxPFP->SetMarkerColor(6);
02112 mHitsUsage = new TH2F("HitsUsage","log10(No.of Used in dE/dx hits) versus log10(Total no. of Tpc Hits",
02113 80,0,8,60,0,6);
02114 }
02115 if (! f && !first) {
02116 for (Int_t i = 0; i < fNZOfBadHits; i++) AddHist(fZOfBadHits[i]);
02117 AddHist(fZOfGoodHits);
02118 AddHist(fPhiOfBadHits);
02119 AddHist(fTracklengthInTpcTotal);
02120 AddHist(fTracklengthInTpc);
02121 AddHist(fTdEdxP70);
02122 AddHist(fTdEdxP70pi);
02123 AddHist(fTdEdxP70e);
02124 AddHist(fTdEdxP70K);
02125 AddHist(fTdEdxP70P);
02126 AddHist(fTdEdxPF);
02127 AddHist(fTdEdxPFpi);
02128 AddHist(fTdEdxPFe);
02129 AddHist(fTdEdxPFK);
02130 AddHist(fTdEdxPFP);
02131 AddHist(mHitsUsage);
02132 }
02133 }
02134 else {
02135 StSPtrVecTrackPidTraits &traits = gTrack->pidTraits();
02136 static StDedxPidTraits *pid, *pid70 = 0, *pidF = 0;
02137 static Double_t TrackLength70, TrackLength, I70, D70, fitZ, fitdZ;
02138 static StProbPidTraits *pidprob = 0;
02139 static Int_t N70, NF;
02140 pid = pid70 = pidF = 0;
02141 TrackLength70 = TrackLength = I70 = D70 = fitZ = fitdZ = 0;
02142 N70 = NF = 0;
02143 for (UInt_t i = 0; i < traits.size(); i++) {
02144 if (! traits[i]) continue;
02145 if ( traits[i]->IsZombie()) continue;
02146 pid = dynamic_cast<StDedxPidTraits*>(traits[i]);
02147 if (pid) {
02148 if (pid->method() == kTruncatedMeanId) {
02149 pid70 = pid; I70 = pid70->mean(); N70 = pid70->numberOfPoints();
02150 TrackLength70 = pid70->length(); D70 = pid70->errorOnMean();
02151 }
02152 if (pid->method() == kLikelihoodFitId) {
02153 pidF = pid;
02154 fitZ = TMath::Log(pidF->mean()+3e-33); NF = pidF->numberOfPoints();
02155 TrackLength = pidF->length(); fitdZ = pidF->errorOnMean();
02156 }
02157 }
02158 else pidprob = dynamic_cast<StProbPidTraits*>(traits[i]);
02159 }
02160 StThreeVectorD g3 = gTrack->geometry()->momentum();
02161 Double_t pMomentum = g3.mag();
02162 if (pid70 && ! pidF) TrackLength = TrackLength70;
02163 if (pid70 && TrackLength70 > 40. && I70 > 0) {
02164 fTdEdxP70->Fill(TMath::Log10(pMomentum), TMath::Log10(I70)+6.);
02165 const StParticleDefinition* pd = gTrack->pidTraits(PidAlgorithm);
02166 if (pd) {
02167 if (TMath::Abs(PidAlgorithm.numberOfSigma(Electron)) < 1) fTdEdxP70e->Fill(TMath::Log10(pMomentum), TMath::Log10(I70)+6.);
02168 if (TMath::Abs(PidAlgorithm.numberOfSigma(Pion)) < 1) fTdEdxP70pi->Fill(TMath::Log10(pMomentum), TMath::Log10(I70)+6.);
02169 if (TMath::Abs(PidAlgorithm.numberOfSigma(Kaon)) < 1) fTdEdxP70K->Fill(TMath::Log10(pMomentum), TMath::Log10(I70)+6.);
02170 if (TMath::Abs(PidAlgorithm.numberOfSigma(Proton)) < 1) fTdEdxP70P->Fill(TMath::Log10(pMomentum), TMath::Log10(I70)+6.);
02171 }
02172 }
02173 if (pidF &&pidprob && TrackLength70 > 40.) {
02174 fTdEdxPF->Fill(TMath::Log10(pMomentum), Log10E*fitZ + 6.);
02175 const StParticleDefinition* pd = gTrack->pidTraits(PidAlgorithm);
02176 if (pd) {
02177 if (pidprob->GetChi2Prob(kPidElectron) > 0.683) fTdEdxPFe->Fill(TMath::Log10(pMomentum), Log10E*fitZ + 6.);
02178 if (pidprob->GetChi2Prob(kPidPion) > 0.683) fTdEdxPFpi->Fill(TMath::Log10(pMomentum), Log10E*fitZ + 6.);
02179 if (pidprob->GetChi2Prob(kPidKaon) > 0.683) fTdEdxPFK->Fill(TMath::Log10(pMomentum), Log10E*fitZ + 6.);
02180 if (pidprob->GetChi2Prob(kPidProton) > 0.683) fTdEdxPFP->Fill(TMath::Log10(pMomentum), Log10E*fitZ + 6.);
02181 }
02182 }
02183 }
02184 first = 2004;
02185 }
02186
02187 void StdEdxY2Maker::BadHit(Int_t iFlag, const StThreeVectorF &xyz) {
02188 if (iFlag >= 0 && iFlag < fNZOfBadHits && fZOfBadHits[iFlag]) fZOfBadHits[iFlag]->Fill(xyz.z());
02189 if (fZOfBadHits[fNZOfBadHits-1]) fZOfBadHits[fNZOfBadHits-1]->Fill(xyz.z());
02190 if (fPhiOfBadHits!= 0) fPhiOfBadHits->Fill(TMath::ATan2(xyz.y(),xyz.x()));
02191 }
02192
02193 void StdEdxY2Maker::Correlations() {
02194
02195 static TH2F *corrI = 0, *corrO = 0, *corrI2 = 0, *corrO2 = 0, *corrI5 = 0, *corrO5 = 0;
02196 static TH2F *corrIw = 0, *corrOw = 0, *corrI2w = 0, *corrO2w = 0, *corrI5w = 0, *corrO5w = 0;
02197 static TH1F *corrI1w = 0, *corrO1w = 0;
02198 if (! corrI) {
02199 gMessMgr->Warning() << "StdEdxY2Maker::Histogramming make Correlation histograms" << endm;
02200 corrI = new TH2F("corrI","Correlation for Inner Sector for pair of nearest rows",
02201 100,-10.,10., 100,-10.,10.);
02202 corrO = new TH2F("corrO","Correlation for Outer Sector for pair of nearest rows",
02203 100,-10.,10., 100,-10.,10.);
02204 corrI2 = new TH2F("corrI2","Correlation for Inner Sector for pair rows & row + 2",
02205 100,-10.,10., 100,-10.,10.);
02206 corrO2 = new TH2F("corrO2","Correlation for Outer Sector for pair rows & row + 2",
02207 100,-10.,10., 100,-10.,10.);
02208 corrI5 = new TH2F("corrI5","Correlation for Inner Sector for pair rows & row + 5",
02209 100,-10.,10., 100,-10.,10.);
02210 corrO5 = new TH2F("corrO5","Correlation for Outer Sector for pair rows & row + 5",
02211 100,-10.,10., 100,-10.,10.);
02212 corrIw = new TH2F("corrIw","Weighted correlation for Inner Sector for pair of nearest rows",
02213 100,-10.,10., 100,-10.,10.);
02214 corrOw = new TH2F("corrOw","Weighted correlation for Outer Sector for pair of nearest rows",
02215 100,-10.,10., 100,-10.,10.);
02216 corrI1w = new TH1F("corrI1w","Weighted distribution for Inner Sector",100,-10.,10.);
02217 corrO1w = new TH1F("corrO1w","Weighted distribution for Outer Sector",100,-10.,10.);
02218 corrI2w = new TH2F("corrI2w","Weighted correlation for Inner Sector for pair rows & row + 2",
02219 100,-10.,10., 100,-10.,10.);
02220 corrO2w = new TH2F("corrO2w","Weighted correlation for Outer Sector for pair rows & row + 2",
02221 100,-10.,10., 100,-10.,10.);
02222 corrI5w = new TH2F("corrI5w","Weighted correlation for Inner Sector for pair rows & row + 5",
02223 100,-10.,10., 100,-10.,10.);
02224 corrO5w = new TH2F("corrO5w","Weighted correlation for Outer Sector for pair rows & row + 5",
02225 100,-10.,10., 100,-10.,10.);
02226 }
02227 for (Int_t k = 0; k < NdEdx; k++) {
02228 Double_t zk = FdEdx[k].zdev;
02229 if (FdEdx[k].Prob > 1.e-12) {
02230 if (FdEdx[k].row > 13) corrO1w->Fill(zk,1./FdEdx[k].Prob);
02231 else corrI1w->Fill(zk,1./FdEdx[k].Prob);
02232 }
02233 for (Int_t m = 0; m < NdEdx; m++){
02234 if (k == m) continue;
02235 Double_t zl = FdEdx[m].zdev;
02236 if (FdEdx[m].row%2 == 1 && FdEdx[m].row - FdEdx[k].row == 1) {
02237 if (FdEdx[k].row > 13) {
02238 corrO->Fill(zk,zl);
02239 if (FdEdx[k].Prob*FdEdx[m].Prob > 1.e-12)
02240 corrOw->Fill(zk,zl,1./(FdEdx[k].Prob*FdEdx[m].Prob));
02241 }
02242 else {
02243 corrI->Fill(zk,zl);
02244 if (FdEdx[k].Prob*FdEdx[m].Prob > 1.e-12)
02245 corrIw->Fill(zk,zl,1./(FdEdx[k].Prob*FdEdx[m].Prob));
02246 }
02247 }
02248 if (FdEdx[m].row%2 == 1 && FdEdx[m].row - FdEdx[k].row == 2) {
02249 if (FdEdx[k].row > 13) {
02250 corrO2->Fill(zk,zl);
02251 if (FdEdx[k].Prob*FdEdx[m].Prob > 1.e-12)
02252 corrO2w->Fill(zk,zl,1./(FdEdx[k].Prob*FdEdx[m].Prob));
02253 }
02254 else {
02255 corrI2->Fill(zk,zl);
02256 if (FdEdx[k].Prob*FdEdx[m].Prob > 1.e-12)
02257 corrI2w->Fill(zk,zl,1./(FdEdx[k].Prob*FdEdx[m].Prob));
02258 }
02259 }
02260 if (FdEdx[m].row%2 == 1 && FdEdx[m].row - FdEdx[k].row == 5) {
02261 if (FdEdx[k].row > 13) {
02262 corrO5->Fill(zk,zl);
02263 if (FdEdx[k].Prob*FdEdx[m].Prob > 1.e-12)
02264 corrO5w->Fill(zk,zl,1./(FdEdx[k].Prob*FdEdx[m].Prob));
02265 }
02266 else {
02267 corrI5->Fill(zk,zl);
02268 if (FdEdx[k].Prob*FdEdx[m].Prob > 1.e-12)
02269 corrI5w->Fill(zk,zl,1./(FdEdx[k].Prob*FdEdx[m].Prob));
02270 }
02271 }
02272 }
02273 }
02274 }
02275
02276 Int_t StdEdxY2Maker::Propagate(const StThreeVectorD &middle,const StThreeVectorD &normal,
02277 const StPhysicalHelixD &helixI, const StPhysicalHelixD &helixO,
02278 Double_t bField,
02279 StThreeVectorD &xyz, StThreeVectorD &dirG, Double_t s[2], Double_t w[2]) {
02280 xyz = StThreeVectorD();
02281 dirG = StThreeVectorD();
02282 s[0] = helixI.pathLength(middle, normal);
02283 s[1] = helixO.pathLength(middle, normal);
02284 Double_t sA[2] = {s[0]*s[0], s[1]*s[1]};
02285 if (sA[0] > 1.e6 || sA[1] > 1.e6) {return 1;}
02286 Double_t sN = sA[0] + sA[1];
02287 w[0] = sA[0]/sN;
02288 w[1] = sA[1]/sN;
02289 if (w[0] > 1.e-4) {xyz += w[0]*helixO.at(s[1]); dirG += w[0]*helixO.momentumAt(s[1],bField);}
02290 if (w[1] > 1.e-4) {xyz += w[1]*helixI.at(s[0]); dirG += w[1]*helixI.momentumAt(s[0],bField);}
02291 return 0;
02292 }
