StRoot: StPass0CalibMaker/StSpaceChargeEbyEMaker.cxx Source File

00001 
00002 //                                                                      //
00003 // StSpaceChargeEbyEMaker performs event-by-event determination         //
00004 // of the space charge correction for tracks, and sets it for           //
00005 // the next event.                                                      //
00006 //                                                                      //
00008 
00009 #include "StSpaceChargeEbyEMaker.h"
00010 #include "StDbUtilities/StMagUtilities.h"
00011 #include "StEventTypes.h"
00012 #include "StMessMgr.h"
00013 #include "StTpcDb/StTpcDb.h"
00014 #include "StTpcDb/StTpcDbMaker.h"
00015 #include "St_db_Maker/St_db_Maker.h"
00016 #include "tables/St_spaceChargeCor_Table.h"
00017 #include "StEvent/StDcaGeometry.h"
00018 #include "StBTofCollection.h"
00019 #include "StBTofHit.h"
00020 #include "StBTofHeader.h"
00021 #include "StBTofPidTraits.h"
00022 #include "StTrackPidTraits.h"
00023 #include "StDetectorDbMaker/St_trigDetSumsC.h"
00024 #include "StEmcUtil/geometry/StEmcGeom.h"
00025 #include "StEmcUtil/projection/StEmcPosition.h"
00026 
00027 #include "TUnixTime.h"
00028 #include "TFile.h"
00029 #include "TH2.h"
00030 #include "TH3.h"
00031 #include "TF1.h"
00032 #include "TNtuple.h"
00033 #include "TPad.h"
00034 
00035 // Histogram ranges:
00036 const int   SCN1 = 200;
00037 const int   SCN2 = 100;
00038 const float SCL = -0.015;
00039 const float SCH =  0.085;
00040 const float SCB = (SCH-SCL)/SCN1;
00041 const float SCX = SCB/TMath::Sqrt(TMath::TwoPi());
00042 
00043 const int   DCN = 125;
00044 const float DCL = -2.5;
00045 const float DCH =  2.5;
00046 
00047 const int   PHN = 72;
00048 const float PI2 = TMath::TwoPi();
00049 
00050 const int   EVN = 1024;
00051 
00052 const int   ZN = 60;
00053 const float ZL = -150.;
00054 const float ZH = 150.;
00055 
00056 const int   GN = 150;
00057 const float GL = -0.3;
00058 const float GH = 1.2;
00059 const int   GZN = 17;
00060 const float GZL = -200.;
00061 const float GZH = 225.;
00062 
00063 const float ZGGRID = 208.707;
00064 
00065 static TF1 ga1("ga1","[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))");
00066 static TF1 ln1("ln1","[0]+((208.707-abs(x))*[1]/100.0)",-150.,150.);
00067 //static TF1 ln1("ln1","[0]+((350.0-abs(x))*[1]/100.0)",-150.,150.);
00068 
00069 ClassImp(StSpaceChargeEbyEMaker)
00070   
00071 //_____________________________________________________________________________
00072 StSpaceChargeEbyEMaker::StSpaceChargeEbyEMaker(const char *name):StMaker(name),
00073     event(0),
00074     Calibmode(kFALSE), PrePassmode(kFALSE), PrePassdone(kFALSE), QAmode(kFALSE),
00075     doNtuple(kFALSE), doReset(kTRUE), doGaps(kFALSE),
00076     inGapRow(0),
00077     vtxEmcMatch(1), vtxTofMatch(0), vtxMinTrks(5),
00078     minTpcHits(25), reqEmcMatch(kTRUE), reqTofMatch(kTRUE),
00079     m_ExB(0),
00080     scehist(0), timehist(0), myhist(0), myhistN(0), myhistP(0),
00081     myhistE(0), myhistW(0), dczhist(0), dcehist(0), dcphist(0),
00082     dcahist(0), dcahistN(0), dcahistP(0), dcahistE(0), dcahistW(0),
00083     gapZhist(0), gapZhistneg(0), gapZhistpos(0), cutshist(0), ntup(0) {
00084 
00085   HN=96;  // max events used, cannot exceed 96 used in header file
00086   MINTRACKS=1500;
00087   //SCALER_ERROR = 0.0006; // by eye from hist: SCvsZDCEpW.gif (liberal)
00088   SCALER_ERROR = 0.0007; // by RMS from hist: SCvsZDCX.gif (liberal)
00089 
00090   // MAXDIFFE is maximum different in sc from last ebye sc
00091   MAXDIFFE =   SCALER_ERROR;
00092   // MAXDIFFA is maximum different in sc from last scaler sc
00093   MAXDIFFA = 2*SCALER_ERROR; // should be about equal to SCALER_ERROR, no?
00094            // Present uncetainties with scalers demands greater tolerance
00095 
00096   runid = 0;
00097   memset(evts,0,HN*sizeof(int));
00098   memset(times,0,HN*sizeof(int));
00099   memset(evtstbin,0,HN*sizeof(float));
00100   evtsnow = 0;
00101 
00102   SetMode(0); // default is mode 0 (no QA, no PrePass)
00103   //DoQAmode(); // For testing
00104 
00105   schist = new TH1F("SpCh","Space Charge",SCN1,SCL,SCH);
00106   schist->SetDirectory(0);
00107   for (int i=0;i<HN;i++){
00108     schists[i] = new TH1F(Form("SpCh%d",i),"Space Charge",SCN1,SCL,SCH);
00109     schists[i]->SetDirectory(0);
00110   }
00111 }
00112 //_____________________________________________________________________________
00113 StSpaceChargeEbyEMaker::~StSpaceChargeEbyEMaker() {
00114   delete schist;
00115   for (int i=0;i<HN;i++) delete schists[i];
00116 }
00117 //_____________________________________________________________________________
00118 Int_t StSpaceChargeEbyEMaker::Finish() {
00119 
00120   if (evt > 0) {
00121     if (PrePassmode) {
00122       if (PrePassdone) WriteTableToFile();
00123       else gMessMgr->Warning("StSpaceChargeEbyEMaker: NO SC => NO OUTPUT");
00124     }
00125     if ((!Calibmode)&&(!PrePassdone)) EvalCalib();
00126     WriteQAHists();
00127   } else {
00128     gMessMgr->Warning("StSpaceChargeEbyEMaker: NO EVENTS => NO OUTPUT");
00129   }
00130 
00131   return kStOk;
00132 }
00133 //_____________________________________________________________________________
00134 Int_t StSpaceChargeEbyEMaker::Init() {
00135 
00136   // Use mode to set switches:
00137   switch (GetMode()) {
00138     case (1) : DoQAmode(); break;
00139     case (2) : DoPrePassmode(); break;
00140     case (3) : DoPrePassmode(); DoQAmode();  break;
00141     case (4) : DoCalib();  break;
00142     case (10): DoNtuple(); break;
00143     case (11): DoNtuple(); DontReset(); break;
00144     case (12): DoNtuple(); DoQAmode(); break;
00145     case (13): DoNtuple(); DontReset(); DoQAmode(); break;
00146     default  : {}
00147   }
00148 
00149   if (Calibmode) doReset = kFALSE;
00150 
00151   evt=0;
00152   oldevt=1;
00153   lastsc=0.;
00154   curhist=0;
00155   lasttime=0;
00156   did_auto=kTRUE;
00157   InitQAHists();
00158   if (QAmode) gMessMgr->Info("StSpaceChargeEbyEMaker: Initializing");
00159   return StMaker::Init();
00160 }
00161 //_____________________________________________________________________________
00162 Int_t StSpaceChargeEbyEMaker::Make() {
00163 
00164   if (QAmode) cutshist->Fill(0);
00165 
00166   // On very first event, determine first event timestamp and
00167   //   set default parameters, unless in Calibmode
00168   if ((!Calibmode) && (tabname.Length() == 0)) SetTableName();
00169   
00170   // Get instance of StMagUtilities
00171   m_ExB = StMagUtilities::Instance();
00172   if (!m_ExB) {
00173     TDataSet *RunLog = GetDataBase("RunLog");
00174     if (!RunLog) gMessMgr->Warning("StSpaceChargeEbyEMaker: No RunLog found.");
00175     m_ExB = new StMagUtilities(gStTpcDb,RunLog,(kSpaceChargeR2 | kGridLeak));
00176   }
00177   lastsc = m_ExB->CurrentSpaceChargeR2();
00178 
00179   // Get StEvent and related info, determine if things are OK
00180   event = (StEvent*) GetInputDS("StEvent");
00181   if (!event) {
00182     gMessMgr->Warning("StSpaceChargeEbyEMaker: no StEvent; skipping event.");
00183     return kStWarn;
00184   }
00185   if (QAmode) cutshist->Fill(1);
00186   // Get runinfo, determine if the magnetic field is nonzero 
00187   // EbyE maker not currently able to handle zero B field
00188   runinfo = event->runInfo();
00189   if ((!runinfo) || (runinfo->magneticField() == 0)) {
00190     gMessMgr->Error("StSpaceChargeEbyEMaker: cannot run due to zero or unknown mag field!");
00191     // Look for any errant zero field SC values as a warning
00192     //   for processing even without EbyE
00193     if ((lastsc != 0) && (runinfo) && (runinfo->magneticField() == 0))
00194       gMessMgr->Warning() << "BE AWARE THAT A NONZERO VALUE OF SPACECHARGE\n"
00195         << "    WAS RETURNED BY DB CALL!\n    (could be a local DB file or"
00196         << " in the actual database)" << endm;
00197     return kStFatal;
00198   }
00199   if (QAmode) cutshist->Fill(2);
00200 
00201   // Select the highest ranked vertex + some quality cuts
00202   StPrimaryVertex* pvtx = event->primaryVertex();
00203   if (!pvtx ||
00204       (pvtx->numberOfDaughters()  < vtxMinTrks)   ||
00205       (pvtx->numMatchesWithBEMC() < vtxEmcMatch)  ||
00206       (pvtx->numMatchesWithBTOF() < vtxTofMatch)) return kStOk;
00207   if (QAmode) cutshist->Fill(3);
00208   StVertexFinderId vid = pvtx->vertexFinderId();
00209   float min_rank = -1e6;
00210   switch (vid) {
00211     case minuitVertexFinder   : min_rank = -5; break;
00212     case ppvVertexFinder      :
00213     case ppvNoCtbVertexFinder : min_rank = 0; break;
00214     default                   : break;
00215   }
00216   if (pvtx->ranking() < min_rank) return kStOk;
00217   if (QAmode) cutshist->Fill(4);
00218   
00219   StSPtrVecTrackNode& theNodes = event->trackNodes();
00220   unsigned int nnodes = theNodes.size();
00221   if (!nnodes) return kStOk;
00222   if (QAmode) cutshist->Fill(5);
00223 
00224   // Store and setup event-wise info
00225   evt++;
00226   int thistime = event->time();
00227   if (lasttime) {
00228     timehist->SetBinContent(evt,thistime-lasttime);
00229   } else {
00230     runid = event->runId();
00231   }
00232   if (doReset) {
00233     if (evt>1) curhist = imodHN(curhist+1);
00234     schists[curhist]->Reset();
00235     if (doGaps) {
00236       gapZhist->Reset();
00237       gapZhistpos->Reset();
00238       gapZhistneg->Reset();
00239     }
00240   } else {
00241     // Do not reset ntuple in calib mode
00242     if (doNtuple && !Calibmode) ntup->Reset();
00243   }
00244 
00245   // Keep time and event number
00246   times[curhist] = thistime;
00247   evts[curhist]=evt;
00248 
00249   // Keep track of # of events in the same time bin
00250   if (thistime == lasttime) evtsnow++;
00251   else evtsnow = 1;
00252   evtstbin[curhist] = evtsnow;
00253 
00254   if (QAmode) {
00255     gMessMgr->Info()
00256       << "used (for this event) SpaceCharge = "
00257       << lastsc << " (" << thistime << ")" << endm;
00258     gMessMgr->Info()
00259       << "zdc west+east = "
00260       << runinfo->zdcWestRate()+runinfo->zdcEastRate() << endm;
00261     gMessMgr->Info()
00262       << "zdc coincidence = "
00263       << runinfo->zdcCoincidenceRate() << endm;
00264   }
00265 
00266   // Fill the StEvent information for the SpaceCharge used in this event
00267   runinfo->setSpaceCharge(lastsc);
00268   runinfo->setSpaceChargeCorrectionMode(m_ExB->GetSpaceChargeMode());
00269   if (!inGapRow) {
00270     if (runinfo->runId() > 10000000) inGapRow = 13; // Run 9+
00271     else if (runinfo->runId() > 0) inGapRow = 12; // Run exists
00272     // else undefined
00273   }
00274  
00275   // Track loop
00276   unsigned int i,j,k;
00277   StThreeVectorD ooo = pvtx->position();
00278 
00279   // Prepare for EMC match
00280   StEmcDetector* bemcDet = 0;
00281   Double_t emcRadius = 0;
00282   static StEmcPosition* emcPosition = 0;
00283   static StEmcGeom* emcGeom = 0;
00284   if (reqEmcMatch) {
00285     bemcDet = event->emcCollection()->detector(kBarrelEmcTowerId);
00286     if (!emcPosition) emcPosition = new StEmcPosition();
00287     if (!emcGeom) emcGeom = StEmcGeom::instance("bemc");
00288     emcRadius = emcGeom->Radius() + 30; // use exit radius, 30cm beyond face
00289   }
00290 
00291   for (i=0; i<nnodes; i++) {
00292       for (j=0; j<theNodes[i]->entries(global); j++) {
00293         if (QAmode) cutshist->Fill(16);
00294         StTrack* tri = theNodes[i]->track(global,j);
00295         if (!tri) continue;
00296         if (QAmode) cutshist->Fill(17);
00297 
00298           const StTrackTopologyMap& map = tri->topologyMap();
00299           //if (! map.trackTpcOnly()) continue;
00300           if (! map.hasHitInDetector(kTpcId)) continue;
00301           if (QAmode) cutshist->Fill(18);
00302           // Multiple silicon hits destroy sDCA <-> SpaceCharge correlation,
00303           // and single hit in SVT is unreliable. Only good config is NO SVT!
00304           if (map.hasHitInDetector(kSvtId)) continue;
00305           if (QAmode) cutshist->Fill(19);
00306           if (map.numberOfHits(kTpcId) < minTpcHits) continue;
00307           if (QAmode) cutshist->Fill(20);
00308 
00309           if (reqTofMatch) {
00310 
00311             const StPtrVecTrackPidTraits& theTofPidTraits = tri->pidTraits(kTofId);
00312             if (!theTofPidTraits.size()) continue;
00313             if (QAmode) cutshist->Fill(21);
00314 
00315             StTrackPidTraits* theSelectedTrait = theTofPidTraits[theTofPidTraits.size()-1];
00316             if (!theSelectedTrait) continue;
00317             if (QAmode) cutshist->Fill(22);
00318             StBTofPidTraits *pidTof = dynamic_cast<StBTofPidTraits *>(theSelectedTrait);
00319             if (!pidTof) continue;
00320             if (QAmode) cutshist->Fill(23);
00321 
00322             int Mflag=pidTof->matchFlag();
00323             //  0: no matching
00324             //  1: 1-1 matching
00325             //  2: 1-2 matching, pick up the one with higher ToT vaule (<25ns) 
00326             //  3: 1-2 matching, pick up the one with closest projection position along local y
00327             if (Mflag <= 0) continue;
00328 
00329           }
00330           if (QAmode) cutshist->Fill(24);
00331 
00332           if (reqEmcMatch) {
00333 
00334             Double_t mEmcThresh = 0.15;
00335             Double_t energyBEMC = -100.0;
00336             UInt_t tower_eta,tower_mod = 0;
00337             Int_t tower_sub = 0;
00338             StThreeVectorD emcTrkMomentum,emcTrkPosition;
00339             if (!(emcPosition->trackOnEmc(&emcTrkPosition,&emcTrkMomentum,
00340                                  tri,runinfo->magneticField()/10.,emcRadius))) continue;
00341             if (QAmode) cutshist->Fill(25);
00342          
00343             Float_t emcEta = emcTrkPosition.pseudoRapidity();
00344             Float_t emcPhi = emcTrkPosition.phi();
00345             Int_t m,e,s,id = 0;
00346             emcGeom->getBin(emcPhi,emcEta,m,e,s);
00347             if (emcGeom->getId(m,e,s,id) == 0) {
00348               tower_mod = m;
00349               tower_eta = e;
00350               tower_sub = s;
00351             }
00352             if (tower_mod < 1 || tower_mod > 120) continue;
00353             if (QAmode) cutshist->Fill(26);
00354  
00355             if (event->emcCollection()) {
00356               StEmcModule* emcMod = bemcDet->module(tower_mod);
00357               StSPtrVecEmcRawHit& emcHits = emcMod->hits();
00358               for (UInt_t emcHit=0; emcHit<emcHits.size(); emcHit++) {
00359                 if ((emcHits[emcHit]) && (emcHits[emcHit]->eta() == tower_eta)
00360                                       && (emcHits[emcHit]->sub() == tower_sub)) {
00361                     energyBEMC = emcHits[emcHit]->energy();
00362                     break; // only one hit
00363                 }
00364               }
00365             }
00366          
00367             if (energyBEMC < mEmcThresh) continue; 
00368 
00369           }
00370           if (QAmode) cutshist->Fill(27);
00371 
00372           StTrackGeometry* triGeom = tri->geometry();
00373 
00374           StThreeVectorF xvec = triGeom->origin();
00375           if (!(xvec.x() || xvec.y() || xvec.z())) continue;
00376           if (QAmode) cutshist->Fill(28);
00377           StThreeVectorF pvec = triGeom->momentum();
00378           if (!(pvec.x() || pvec.y())) continue;
00379           if (QAmode) cutshist->Fill(29);
00380 
00381           float oldPt = pvec.perp();
00382           if (oldPt < 0.0001) continue;
00383           if (QAmode) cutshist->Fill(30);
00384 
00385           int e_or_w = 0; // east is -1, west is +1
00386           if (pvec.z() * xvec.z() > 0) e_or_w = ( (xvec.z() > 0) ? 1 : -1 );
00387 
00388           StPhysicalHelixD hh = triGeom->helix();
00389 
00390           Float_t eta=pvec.pseudoRapidity();
00391           Float_t phi=0;
00392           Double_t pathlen = 0;
00393           //Float_t DCA=hh.geometricSignedDistance(0,0); // for testing only
00394           Float_t DCA3=-999;
00395           Float_t DCA2=-999;
00396           Double_t DCAerr = 0.;
00397           StDcaGeometry* triDcaGeom = ((StGlobalTrack*) tri)->dcaGeometry();
00398           if (triDcaGeom) {
00399             StPhysicalHelixD dcahh = triDcaGeom->helix();
00400             DCA3 = dcahh.distance(ooo,kFALSE);
00401             DCA2 = dcahh.geometricSignedDistance(ooo.x(),ooo.y());
00402             // helix() gets the sign of DCA2, thelix() gets the error
00403             THelixTrack thelix = triDcaGeom->thelix();
00404             thelix.Dca(ooo.x(),ooo.y(),&DCAerr);
00405             phi = TMath::ATan2(dcahh.cy(pathlen),dcahh.cx(pathlen));
00406           } else {
00407             DCA3 = hh.distance(ooo,kFALSE);
00408             DCA2 = hh.geometricSignedDistance(ooo.x(),ooo.y());
00409             pathlen = hh.pathLength(ooo.x(),ooo.y());
00410             phi = TMath::ATan2(hh.cy(pathlen),hh.cx(pathlen));
00411           }
00412           if (DCA3 > 4) continue; // cut out pileup tracks!
00413           if (QAmode) cutshist->Fill(31);
00414           Int_t ch = (int) triGeom->charge();
00415 
00416           Int_t PCT = 0;
00417           Float_t rerrors[64];
00418           Float_t rphierrors[64];
00419           memset(rerrors,64*sizeof(Float_t),0);
00420           memset(rphierrors,64*sizeof(Float_t),0);
00421           StPtrVecHit& hits = tri->detectorInfo()->hits();
00422           for (k=0;k<hits.size();k++) {
00423             StHit* hit = hits[k];
00424             unsigned int maskpos = 0;
00425             switch (hit->detector()) {
00426               case (kTpcId) :
00427                 if ((hit->position().z() > 1 && ((StTpcHit*) hit)->sector() > 12) ||
00428                     (hit->position().z() <-1 && ((StTpcHit*) hit)->sector() < 13)) PCT++;
00429                 maskpos = 7 + ((StTpcHit*) hit)->padrow(); break;
00430               case (kSvtId) :
00431                 maskpos = ((StSvtHit*) hit)->layer(); break;
00432               case (kSsdId) :
00433                 maskpos = 7; break;
00434               default :
00435                 maskpos = 0;
00436             }
00437             if (maskpos) {
00438               StThreeVectorF herrVec = hit->positionError();
00439               Float_t herr = herrVec.perp();
00440               rerrors[maskpos] = herr;
00441               rphierrors[maskpos] = herr;
00442             }
00443           }
00444           if (PCT) continue; // Track has post-crossing hits
00445           if (QAmode) cutshist->Fill(32);
00446           
00447           Float_t space = 10000.;
00448           if (!(m_ExB->PredictSpaceChargeDistortion(ch,oldPt,ooo.z(),
00449              eta,phi,DCA2,map.data(0),map.data(1),rerrors,rphierrors,space))) continue;
00450           if (QAmode) cutshist->Fill(33);
00451 
00452           Double_t spaceErr = TMath::Abs(space*DCAerr/DCA2);
00453           space += lastsc;  // Assumes additive linearity of space charge!
00454           if (spaceErr > 0) {
00455             // Fill SpaceCharge accounting for prediction error
00456             ga1.SetParameters(SCX/spaceErr,space,spaceErr);
00457             for (Int_t si=1;si<=SCN1;si++) {
00458               Double_t sx = schists[curhist]->GetBinCenter(si);
00459               if (TMath::Abs((sx-space)/spaceErr) < 4.5) {
00460                 // only within +/-4.5 sigma
00461                 schists[curhist]->Fill(sx,ga1.Eval(sx));
00462               }
00463             }
00464           } else {
00465             schists[curhist]->Fill(space);
00466           }
00467           FillQAHists(DCA2,space,ch,hh,e_or_w);
00468 
00469 
00470           if ((doGaps) &&
00471               (map.hasHitInRow(kTpcId,inGapRow))&&(map.hasHitInRow(kTpcId,14)) &&
00472               (map.hasHitInRow(kTpcId,11))&&(map.hasHitInRow(kTpcId,15)) &&
00473               (e_or_w!=0) && (TMath::Abs(ch)==1) && (oldPt>0.3))
00474             FillGapHists(tri,hh,e_or_w,ch);
00475 
00476       } // loop over j tracks
00477   } // loop over i Nodes
00478   if (QAmode) cutshist->Fill(6);
00479 
00480 
00481   ntrks[curhist] = schists[curhist]->Integral();
00482 
00483   // Wrap it up and make a decision
00484   int result = DecideSpaceCharge(thistime);
00485 
00486   if (doGaps) DetermineGaps();
00487   if (doNtuple) {
00488       static float X[39];
00489       static float ntent = 0.0;
00490       static float nttrk = 0.0;
00491 
00492       if (ntent == 0.0) for (i=0; i<39; i++) X[i] = 0.0;
00493       ntent++;  // # entries since last reset, including this one
00494       float last_nttrk = nttrk;
00495       nttrk = ntrks[curhist];  // # tracks since last reset, including these
00496       float s0 = ( nttrk ? last_nttrk / nttrk : 0 );
00497       float s1 = 1.0 - s0; // fraction of tracks from current event
00498 
00499       if (QAmode) {
00500         gMessMgr->Info() << "reset = " << doReset << endm;
00501         gMessMgr->Info() << "nevts = " << ntent << endm;
00502         gMessMgr->Info() << "ntrks = " << nttrk << endm;
00503       }
00504 
00505       float ee;
00506       int fbin = evt + 1 - ((int) ntent);
00507 
00508       X[0]  = sc;
00509       X[1]  = FindPeak(dcehist->ProjectionY("_dy",fbin,evt),ee);
00510       X[2]  = s0*X[2] + s1*runinfo->zdcCoincidenceRate();
00511       X[3]  = s0*X[3] + s1*runinfo->zdcWestRate();
00512       X[4]  = s0*X[4] + s1*runinfo->zdcEastRate();
00513       X[5]  = s0*X[5] + s1*runinfo->bbcCoincidenceRate();
00514       X[6]  = s0*X[6] + s1*runinfo->bbcWestRate();
00515       X[7]  = s0*X[7] + s1*runinfo->bbcEastRate();
00516       X[8]  = s0*X[8] + s1*runinfo->bbcBlueBackgroundRate();
00517       X[9]  = s0*X[9] + s1*runinfo->bbcYellowBackgroundRate();
00518       X[10] = s0*X[10] + s1*runinfo->initialBeamIntensity(blue);   // west-bound
00519       X[11] = s0*X[11] + s1*runinfo->initialBeamIntensity(yellow); // east-bound
00520       X[12] = runinfo->beamFillNumber(blue);
00521       X[13] = runinfo->magneticField();
00522       X[14] = event->runId();
00523       X[15] = event->id();
00524       //X[16] = sceast;
00525       //X[17] = scwest;
00526       if ((QAmode) && (evt <= EVN)) {
00527         X[16] = FindPeak(dcahistN->ProjectionZ("_dnz",fbin,evt,1,PHN),ee);
00528         X[17] = FindPeak(dcahistP->ProjectionZ("_dpz",fbin,evt,1,PHN),ee);
00529         X[18] = FindPeak(dcahistE->ProjectionZ("_dez",fbin,evt,1,PHN),ee);
00530         X[19] = FindPeak(dcahistW->ProjectionZ("_dwz",fbin,evt,1,PHN),ee);
00531       }
00532       X[20] = gapZfitslope;
00533       X[21] = gapZfitintercept;
00534       X[22] = gapZdivslope;
00535       X[23] = gapZfitslopeneg;
00536       X[24] = gapZfitinterceptneg;
00537       X[25] = gapZdivslopeneg;
00538       X[26] = gapZfitslopepos;
00539       X[27] = gapZfitinterceptpos;
00540       X[28] = gapZdivslopepos;
00541       X[29] = gapZfitslopeeast;
00542       X[30] = gapZfitintercepteast;
00543       X[31] = gapZdivslopeeast;
00544       X[32] = gapZfitslopewest;
00545       X[33] = gapZfitinterceptwest;
00546       X[34] = gapZdivslopewest;
00547       X[35] = s0*X[35] + s1*runinfo->spaceCharge();
00548       X[36] = s0*X[36] + s1*((float) (runinfo->spaceChargeCorrectionMode()));
00549       X[37] = s0*X[37] + s1*St_trigDetSumsC::Nc(runinfo->zdcCoincidenceRate(),
00550                                  runinfo->zdcEastRate(),runinfo->zdcWestRate());
00551       X[38] = s0*X[38] + s1*St_trigDetSumsC::Nc(runinfo->bbcCoincidenceRate(),
00552                                  runinfo->bbcEastRate(),runinfo->bbcWestRate());
00553 
00554       // In calib mode, only fill when doReset (we found an sc)
00555       if (doReset || !Calibmode) ntup->Fill(X);
00556 
00557       if (doReset) {ntent = 0.0; nttrk = 0.0; }
00558 
00559   }
00560           
00561   lasttime = thistime;
00562 
00563   return result;
00564 }
00565 //_____________________________________________________________________________
00566 Int_t StSpaceChargeEbyEMaker::DecideSpaceCharge(int time) {
00567 
00568   // curhist has only this event
00569   // curhist-1 has past two events...
00570   // curhist-x has past (x+1) events
00571   // curhist-(HN-1) == curhist+1 has past HN events
00572 
00573   Bool_t QAout = QAmode || PrePassmode;
00574   Bool_t do_auto = kTRUE;
00575   Bool_t few_stats = kTRUE;
00576   Bool_t large_err = kTRUE;
00577   Bool_t large_dif = kTRUE;
00578 
00579   // Cuts on difference from previous sc measure:
00580   // If last event was auto, use MAXDIFFA, else use between MAXDIFFE & MAXDIFFA
00581   //   scaled by oldness of previous sc measure (curhist-1)
00582   float maxdiff,dif=0;
00583   if (did_auto)
00584     maxdiff = MAXDIFFA;
00585   else
00586     maxdiff = MAXDIFFA - (MAXDIFFA-MAXDIFFE)*oldness(imodHN(curhist-1));
00587 
00588   // More than 30 seconds since last used event? Forget it...
00589   int timedif = time-lasttime;
00590   if (QAout) {
00591     gMessMgr->Info() << "time since last event = "
00592       << timedif << endm;
00593     gMessMgr->Info() << "curhist = "
00594       << curhist << endm;
00595   }
00596   float ntrkstot = 0; // running sum using oldness scale factor
00597   Bool_t decideFromData = ((PrePassmode) || (Calibmode) || (lasttime==0) || (timedif < 30));
00598   if (decideFromData) {
00599   
00600     int isc;
00601     static int iscMax = 1;  // use only one hist for calib mode, and...
00602     if (!Calibmode && iscMax<HN) iscMax = curhist+1; // don't use uninitialized
00603     for (isc=0; isc<iscMax; isc++) {
00604       int i = imodHN(curhist-isc);
00605       ntrkstot += ntrks[i] * oldness(i);
00606       if (QAout) {
00607         if (!isc) gMessMgr->Info("Building with: i, ni, oi, nt:");
00608         gMessMgr->Info() << "Building with: " << i << ", "
00609           << ntrks[i] << ", " << oldness(i) << ", " << ntrkstot << endm;
00610       }
00611 
00612       // Too little data collected? Keep trying...
00613       few_stats = ntrkstot < MINTRACKS;
00614       if (!few_stats) {
00615         BuildHist(i);
00616         FindSpaceCharge();
00617         if (QAout) gMessMgr->Info()
00618           << "sc = " << sc << " +/- " << esc << endm;
00619         large_err = (esc == 0) || (esc > SCALER_ERROR);
00620         if (!large_err) {
00621           if (PrePassmode) { do_auto=kFALSE; break; }
00622           // otherwise, check for big jumps
00623           dif = TMath::Abs(sc-lastsc);
00624           large_dif = dif > maxdiff;
00625           if (!large_dif || Calibmode) {
00626             oldevt = evts[i];
00627             do_auto=kFALSE;
00628             break;
00629           }
00630         }
00631       }
00632 
00633       // shouldn't need to go past oldest event previously used?
00634       // tough to know - allow for now as of 11 Jan 2008
00635       // if (evts[i] <= oldevt)  break;
00636     }
00637     if (QAout && (isc == HN)) gMessMgr->Info()
00638       << "STORED DATA EXHAUSTED: "
00639       << HN << " events" << endm;
00640   }
00641 
00642   did_auto = do_auto;
00643 
00644   // In normal  mode, do_auto decides whether to use automatic SC from DB
00645   // In PrePass mode, do_auto decides when we're ready to stop
00646   // In Calib   mode, do_auto decides when to save entries and reset
00647 
00648   if (do_auto) {
00649     if (QAout && decideFromData) {
00650       if (few_stats) gMessMgr->Info()
00651         << "(RECENT) STATS TOO FEW: "
00652         << ntrkstot << " (" << MINTRACKS << ")" << endm;
00653       else if (large_err) gMessMgr->Info()
00654         << "FIT ERROR TOO LARGE: "
00655         << esc << " (" << SCALER_ERROR << ")" << endm;
00656       else if (large_dif) gMessMgr->Info()
00657         << "DIFFERENCE TOO LARGE: "
00658         << dif << " (" << maxdiff << ")" << endm;
00659     }
00660     gMessMgr->Info("using auto SpaceCharge");
00661     if (Calibmode) doReset = kFALSE;
00662     else m_ExB->AutoSpaceChargeR2();
00663   } else {
00664     gMessMgr->Info() << "using SpaceCharge = "
00665       << sc << " +/- " << esc << " (" << time << ")" << endm;
00666     scehist->SetBinContent(evt,sc);
00667     scehist->SetBinError(evt,esc);
00668     if (PrePassmode) {
00669       PrePassdone = kTRUE;
00670       return kStStop; // We're happy! Let's stop!
00671     }
00672     if (Calibmode) {
00673       doReset = kTRUE;
00674       return kStStop; // We're happy? Let's stop!
00675     }
00676     else m_ExB->ManualSpaceChargeR2(sc,m_ExB->CurrentSpaceChargeEWRatio());
00677   }
00678   return kStOk;
00679 }
00680 //_____________________________________________________________________________
00681 void StSpaceChargeEbyEMaker::FindSpaceCharge() {
00682   esc = 0.;
00683   double res = FindPeak(schist,esc);
00684   if (res > -500.) sc = res;
00685 }
00686 //_____________________________________________________________________________
00687 double StSpaceChargeEbyEMaker::FindPeak(TH1* hist,float& pkwidth) {
00688 
00689   pkwidth = 0.;
00690   if (hist->Integral() < 100.0) return -998.;
00691   double mn = hist->GetMean();
00692   double rms = TMath::Abs(hist->GetRMS());
00693   mn *= 1.1; rms *= 1.5;
00694   double lr = mn-rms;
00695   double ur = mn+rms;
00696   double pmax = TMath::Max(hist->GetMaximum(),0.);
00697   double lp = pmax*0.001;
00698   double up = pmax*10.0;
00699   ga1.SetParameters(pmax,mn,rms*0.5);
00700   ga1.SetRange(lr,ur);
00701   ga1.SetParLimits(0,lp,up); // Loglikelihood only works with positive functions
00702   int fitResult = hist->Fit(&ga1,
00703     (gROOT->GetVersionInt() >= 53000 ? "WLR0Q" : "LLR0Q")); // Loglikelihood options changed!
00704   if (fitResult != 0) return -999.;
00705   double rp = ga1.GetParameter(0);
00706   if (rp == lp || rp == up) return -998;
00707   pkwidth = TMath::Abs(ga1.GetParError(1));
00708   return ga1.GetParameter(1);
00709 
00710 }
00711 //_____________________________________________________________________________
00712 void StSpaceChargeEbyEMaker::InitQAHists() {
00713 
00714   scehist  = new TH1F("SpcChgEvt","SpaceCharge fit vs. Event",
00715                       EVN,0.,EVN);
00716   timehist = new TH1F("EvtTime","Event Times",
00717                       EVN,0.,EVN);
00718   dcehist  = new TH2F("DcaEve","psDCA vs. Event",
00719                       EVN,0.,EVN,DCN,DCL,DCH);
00720   dcphist  = new TH3F("DcaPhi","psDCA vs. Phi",
00721                       PHN,0,PI2,DCN,DCL,DCH,(QAmode ? 3 : 1),-1.5,1.5);
00722 
00723   AddHist(scehist);
00724   AddHist(timehist);
00725   AddHist(dcehist);
00726   AddHist(dcphist);
00727 
00728   if (QAmode) {
00729     myhist   = new TH3F("SpcEvt","SpaceCharge vs. Phi vs. Event",
00730                         EVN,0.,EVN,PHN,0,PI2,SCN2,SCL,SCH);
00731     dcahist  = new TH3F("DcaEvt","psDCA vs. Phi vs. Event",
00732                         EVN,0.,EVN,PHN,0,PI2,DCN,DCL,DCH);
00733     dczhist  = new TH2F("DcaZ","psDCA vs. Z",
00734                         //80,-200,200,250,-5.0,5.0);
00735                         ZN,ZL,ZH,DCN,DCL,DCH);
00736     myhistN  = new TH3F("SpcEvtN","SpaceCharge vs. Phi vs. Event Neg",
00737                         EVN,0.,EVN,PHN,0,PI2,SCN2,SCL,SCH);
00738     myhistP  = new TH3F("SpcEvtP","SpaceCharge vs. Phi vs. Event Pos",
00739                         EVN,0.,EVN,PHN,0,PI2,SCN2,SCL,SCH);
00740     myhistE  = new TH3F("SpcEvtE","SpaceCharge vs. Phi vs. Event East",
00741                         EVN,0.,EVN,PHN,0,PI2,SCN2,SCL,SCH);
00742     myhistW  = new TH3F("SpcEvtW","SpaceCharge vs. Phi vs. Event West",
00743                         EVN,0.,EVN,PHN,0,PI2,SCN2,SCL,SCH);
00744     dcahistN = new TH3F("DcaEvtN","psDCA vs. Phi vs. Event Neg",
00745                         EVN,0.,EVN,PHN,0,PI2,DCN,DCL,DCH);
00746     dcahistP = new TH3F("DcaEvtP","psDCA vs. Phi vs. Event Pos",
00747                         EVN,0.,EVN,PHN,0,PI2,DCN,DCL,DCH);
00748     dcahistE = new TH3F("DcaEvtE","psDCA vs. Phi vs. Event East",
00749                         EVN,0.,EVN,PHN,0,PI2,DCN,DCL,DCH);
00750     dcahistW = new TH3F("DcaEvtW","psDCA vs. Phi vs. Event West",
00751                         EVN,0.,EVN,PHN,0,PI2,DCN,DCL,DCH);
00752     cutshist = new TH1I("CutsHist","Step at which tracks and events are cut",
00753                         64,-0.5,63.5); // Use < 16 for event-wise cuts
00754     AddHist(myhist);
00755     AddHist(dcahist);
00756     AddHist(dczhist);
00757     AddHist(myhistN);
00758     AddHist(myhistP);
00759     AddHist(myhistE);
00760     AddHist(myhistW);
00761     AddHist(dcahistN);
00762     AddHist(dcahistP);
00763     AddHist(dcahistE);
00764     AddHist(dcahistW);
00765     AddHist(cutshist);
00766   }
00767 
00768   if (doNtuple) ntup = new TNtuple("SC","Space Charge",
00769     "sc:dca:zdcx:zdcw:zdce:bbcx:bbcw:bbce:bbcbb:bbcyb:intb:inty:fill:mag:run:event:dcan:dcap:dcae:dcaw:gapf:gapi:gapd:gapfn:gapin:gapdn:gapfp:gapip:gapdp:gapfe:gapie:gapde:gapfw:gapiw:gapdw:usc:uscmode:zdcc:bbcc");
00770     //"sc:dca:zdcx:zdcw:zdce:bbcx:bbcw:bbce:bbcbb:bbcyb:intb:inty:fill:mag:run:event:sce:scw:dcae:dcaw:gapf:gapd:gapfn:gapdn:gapfp:gapdp:gapfe:gapde:gapfw:gapdw:usc:uscmode");
00771 
00772   if (doGaps) {
00773     gapZhist = new TH2F("Gaps","Gaps",GN,GL,GH,GZN,GZL,GZH);
00774     gapZhistneg = new TH2F("Gapsneg","Gaps Neg",GN,GL,GH,GZN,GZL,GZH);
00775     gapZhistpos = new TH2F("Gapspos","Gaps Pos",GN,GL,GH,GZN,GZL,GZH);
00776   }
00777 
00778 }
00779 //_____________________________________________________________________________
00780 void StSpaceChargeEbyEMaker::WriteQAHists() {
00781 // Only if QAmode or doNtuple or doGaps
00782 
00783   if (!(QAmode || doNtuple || doGaps)) return;
00784 
00785   if (runid == 0) {
00786     gMessMgr->Error("StSpaceChargeEbyEMaker: No runid => No output");
00787     return;
00788   }
00789 
00790   const char* f1 = GetName();
00791   runid -= 1000000*(runid/1000000); // Remove the year, optional
00792   runid *= 100; // limits one run to 100 files!
00793   TString fname = "./hists";
00794   if (PrePassmode) fname.Append("Pre");
00795   gSystem->cd(fname.Data());
00796   while (f1) {
00797     fname = Form("Hist%d.root",runid);
00798     f1 = gSystem->Which(".",fname.Data());
00799     runid++;
00800   }
00801 
00802   TFile ff(fname.Data(),"RECREATE");
00803   if (QAmode) {
00804     myhist->Write();
00805     dcehist->Write();
00806     dcphist->Write();
00807     dcahist->Write();
00808     dczhist->Write();
00809     myhistN->Write();
00810     dcahistN->Write();
00811     myhistP->Write();
00812     dcahistP->Write();
00813     myhistE->Write();
00814     dcahistE->Write();
00815     myhistW->Write();
00816     dcahistW->Write();
00817     scehist->Write();
00818     timehist->Write();
00819     cutshist->Write();
00820   }
00821   if (doGaps) {
00822     gapZhist->Write();
00823     gapZhistneg->Write();
00824     gapZhistpos->Write();
00825   }
00826   if (doNtuple) ntup->Write();
00827   ff.Close();
00828 
00829   gMessMgr->Info() << "QA hists file: " << fname.Data() << endm;
00830 
00831   gSystem->cd("..");
00832 
00833 }
00834 //_____________________________________________________________________________
00835 void StSpaceChargeEbyEMaker::FillQAHists(float DCA, float space, int ch,
00836                                          StPhysicalHelixD& hh, int e_or_w) {
00837   // Find a "Phi" for the track
00838   pairD pls = hh.pathLength(97.0);
00839   double pl = pls.first;
00840   if (TMath::Abs(pls.second) < TMath::Abs(pls.first)) pl = pls.second;
00841   StThreeVectorD hh_at_pl = hh.at(pl);
00842   float Phi = hh_at_pl.phi();
00843   while (Phi < 0) Phi += PI2;
00844   while (Phi >= TMath::TwoPi()) Phi -= PI2;
00845   
00846   // To pile all sectors atop each other:
00847   // while (Phi >= TMath::Pi()/6.) Phi -= TMath::Pi()/6.;
00848 
00849   float evtn = ((float) evt) - 0.5;
00850   dcehist->Fill(evtn,DCA);
00851   dcphist->Fill(Phi,DCA,(float) e_or_w);
00852 
00853   if (QAmode) {
00854     myhist->Fill(evtn,Phi,space);
00855     dcahist->Fill(evtn,Phi,DCA);
00856     if (ch > 0) {
00857       myhistP->Fill(evtn,Phi,space);
00858       dcahistP->Fill(evtn,Phi,DCA);
00859     } else {
00860       myhistN->Fill(evtn,Phi,space);
00861       dcahistN->Fill(evtn,Phi,DCA);
00862     }
00863     if (e_or_w > 0) {
00864       myhistW->Fill(evtn,Phi,space);
00865       dcahistW->Fill(evtn,Phi,DCA);
00866     } else if (e_or_w < 0) {
00867       myhistE->Fill(evtn,Phi,space);
00868       dcahistE->Fill(evtn,Phi,DCA);
00869     }
00870     if ((e_or_w != 0) && (TMath::Abs(hh.dipAngle()) < 0.05)) dczhist->Fill(hh_at_pl.z(),DCA);
00871   }
00872 }
00873 //_____________________________________________________________________________
00874 int StSpaceChargeEbyEMaker::imodHN(int i) {
00875   // Keep index in bounds of circular queue
00876   return ( i >= HN ? imodHN(i-HN) : (i < 0 ? imodHN(i+HN) : i) );
00877 }
00878 //_____________________________________________________________________________
00879 float StSpaceChargeEbyEMaker::oldness(int i, int j) {
00880   // Deterime how to treat relative "age" of event
00881   // In PrePassmode, earliest events are most important!
00882   static float decay_const = -0.12;
00883   //static float decay_const = -0.15;
00884   float s = 1.0;
00885   if (!PrePassmode) { // Weight newest the most (or evenly for PrePass)
00886     if (j<0) j = curhist;
00887 
00888     // Weight in sub-second intervals by # of events because
00889     // times have only 1 second granularity (best we can do).
00890     // Method assumes time-ordering of events, which is violated
00891     // perhaps only occasionally from DAQ.
00892     int k = i;
00893     while (k!=j) {
00894       k = imodHN(k+1);
00895       if (times[k] != times[i]) { k = imodHN(k-1); break; }
00896     }
00897     // # seconds + fraction of a second:
00898     float time_factor = (times[j]-times[i]) + (1.-(evtstbin[i]/evtstbin[k]));
00899     //float time_factor = (times[j]-times[i]);
00900     s = exp( decay_const * time_factor );
00901   }
00902   return s;
00903 }
00904 //_____________________________________________________________________________
00905 void StSpaceChargeEbyEMaker::BuildHist(int i) {
00906   // Build up one histogram from several events
00907   schist->Reset();
00908   int isc = curhist;
00909   schist->Add(schists[isc],1.0);
00910   while (isc != i) {
00911     isc = imodHN(isc-1);
00912     schist->Add(schists[isc],oldness(isc));
00913   }
00914 }
00915 //_____________________________________________________________________________
00916 void StSpaceChargeEbyEMaker::SetTableName() {
00917   // Problem caused if first event comes later in time than other events.
00918   // Solution: subtract 10 seconds...
00919   TUnixTime ux(GetDateTime(),1);
00920   ux+=-10;
00921   int date,time;
00922   ux.GetGTime(date,time);
00923   gMessMgr->Info() << "first event date = " << date << endm;
00924   gMessMgr->Info() << "first event time = " << time << endm;
00925   tabname = Form("./StarDb/Calibrations/rich/spaceChargeCorR2.%08d.%06d.C",date,time);
00926 
00927   // Set default parameters based on data time
00928   //   for non-calib modes
00929   if (date < 20071000) {
00930     setVtxEmcMatch(0);
00931     setReqEmcMatch(kFALSE);
00932     setVtxTofMatch(0);
00933     setReqTofMatch(kFALSE);
00934     setVtxMinTrks(10);
00935   } else if (date < 20090000) {
00936     setVtxEmcMatch(1);
00937     setReqEmcMatch(kFALSE);
00938     setVtxTofMatch(0);
00939     setReqTofMatch(kFALSE);
00940     setVtxMinTrks(5);
00941   }
00942 }
00943 //_____________________________________________________________________________
00944 void StSpaceChargeEbyEMaker::WriteTableToFile(){
00945   gMessMgr->Info() << "Writing new table to:\n  "
00946     << tabname.Data() << endm;
00947   TString dirname = gSystem->DirName(tabname.Data());
00948   TString estr = dirname;
00949   estr.Prepend("mkdir -p ");
00950   gSystem->Exec(estr.Data());
00951   if (gSystem->OpenDirectory(dirname.Data())==0) {
00952     if (gSystem->mkdir(dirname.Data())) {
00953       gMessMgr->Warning() << "Directory creation failed for:\n  " << dirname
00954         << "\n  Putting table file in current directory" << endm;
00955       tabname.Remove(0,tabname.Last('/')).Prepend(".");
00956     }
00957   }
00958   ofstream *out = new ofstream(tabname.Data());
00959   SCTable()->SavePrimitive(*out,"");
00960   return;
00961 }
00962 //_____________________________________________________________________________
00963 St_spaceChargeCor* StSpaceChargeEbyEMaker::SCTable() {
00964   St_spaceChargeCor* table = new St_spaceChargeCor("spaceChargeCorR2",1);
00965   spaceChargeCor_st* row = table->GetTable();
00966   row->fullFieldB = sc;
00967   row->halfFieldB = sc;
00968   row->zeroField  = (float) evt;
00969   row->halfFieldA = sc;
00970   row->fullFieldA = sc;
00971   row->satRate = 1.0;
00972   row->factor = 1.0;
00973   row->detector = 3;
00974   row->offset = 0;
00975   row->ewratio = m_ExB->CurrentSpaceChargeEWRatio();
00976   table->SetNRows(1);
00977   return table;
00978 }
00979 //_____________________________________________________________________________
00980 float StSpaceChargeEbyEMaker::FakeAutoSpaceCharge() {
00981   // Use this to mimic what is done with the scalers, using your own code
00982   float zdcsum = runinfo->zdcWestRate()+runinfo->zdcEastRate();
00983   float sc = 6e-8 * zdcsum;
00984   
00985   return sc;
00986 }
00987 //_____________________________________________________________________________
00988 void StSpaceChargeEbyEMaker::FillGapHists(StTrack* tri, StPhysicalHelixD& hh,
00989                                           int e_or_w, int ch) {
00990   float fsign = ( event->runInfo()->magneticField() < 0 ? -1 : 1 );
00991   StPtrVecHit hts = tri->detectorInfo()->hits(kTpcId);
00992   float gap = 0.; float zgap = 0.; int ct=0;
00993   for (UInt_t ht=0; ht<hts.size(); ht++) {
00994     StTpcHit* hit = (StTpcHit*) hts[ht];
00995     UInt_t prow = hit->padrow();
00996     if ((prow != inGapRow) && (prow != 14)) continue;
00997     float gsign = ( prow == 14 ? -1 : 1 );
00998     const StThreeVectorF& hp = hit->position();
00999 
01000     // Avoid sector edges
01001     float hphi = hp.phi() + TMath::TwoPi();
01002     while (hphi > TMath::Pi()/12.) hphi -= TMath::Pi()/6.;
01003     if (TMath::Abs(hphi) > 0.75*TMath::Pi()/12.) break;
01004 
01005     if (inGapRow==13) zgap += (hp.z() / 7.595) * ( prow == 14 ? 2.2 : 5.395 ); // ~z at gap
01006     else if (inGapRow==12) zgap += (hp.z() / 12.795) * ( prow == 14 ? 7.4 : 5.395 ); // ~z at gap
01007     else return;
01008 
01009     // Measurement method described at:
01010     // http://drupal.star.bnl.gov/STAR/blog/genevb/2010/feb/21/gridleak-update-using-residuals-along-padrows
01011     Double_t residual = hh.geometricSignedDistance(hp.x(),hp.y());
01012 
01013     Int_t sec = hit->sector();
01014     Double_t sector_angle = (TMath::Pi()/6.) * (sec < 13 ? 3 - sec : sec - 21);
01015     Double_t pathlen = hh.pathLength(hp.x(),hp.y());
01016     Double_t theta = TMath::ATan2(hh.cy(pathlen),hh.cx(pathlen)) - sector_angle;
01017     Double_t phi = hp.phi() - sector_angle;
01018 
01019     Double_t Eff1 = TMath::Cos(theta);
01020     Double_t Eff2 = 1;
01021     Double_t Eff3 = 0;
01022 
01023     Float_t x1[3],x2[3];
01024     x1[0] = hp.perp()*TMath::Sin(phi);
01025     x1[1] = hp.perp()*TMath::Cos(phi);
01026     x1[2] = hp.z();
01027     m_ExB->Undo3DGridLeakDistortion(x1,x2,sec);
01028     Double_t dX = x2[0]-x1[0];
01029     if (TMath::Abs(dX) > 1e-20) {
01030       // warning: no available GridLeak calculation may lead
01031       //   to slightly different results
01032       Eff3 = gsign * ((x2[1]-x1[1])/dX) * TMath::Sin(theta);
01033       x1[0] = 0;
01034       m_ExB->Undo3DGridLeakDistortion(x1,x2,sec);
01035       Eff2 = (x2[0]-x1[0])/dX;
01036     }
01037 
01038     Double_t DistortionX = Eff2 * residual / (Eff1 + Eff3);
01039  
01040     gap += fsign * gsign * DistortionX;
01041     ct++;
01042   }
01043 
01044   float abs_zgap = TMath::Abs(zgap);
01045   if ((ct==2) && (abs_zgap<200.0) && (abs_zgap>10.0)) {
01046      gapZhist->Fill(gap,zgap);
01047      if (ch==1) gapZhistpos->Fill(gap,zgap);
01048      else gapZhistneg->Fill(gap,zgap);
01049      if (abs_zgap<150 && abs_zgap>25) { // Restrict the Z range further
01050        // normalize at z=100cm from ggrid
01051        float gap_scaled = (gap * 100.0) / (ZGGRID - abs_zgap);
01052        //float gap_scaled = (gap * 100.0) / (350.0 - abs_zgap);
01053        float z_beyond = ZGGRID+1.0;
01054        gapZhist->Fill(gap_scaled,z_beyond);
01055        if (ch==1) gapZhistpos->Fill(gap_scaled,z_beyond);
01056        else gapZhistneg->Fill(gap_scaled,z_beyond);
01057      }
01058    }
01059 }
01060 //_____________________________________________________________________________
01061 void StSpaceChargeEbyEMaker::DetermineGaps() {
01062   DetermineGapHelper(gapZhistneg,gapZfitslopeneg,gapZfitinterceptneg,gapZdivslopeneg);
01063   DetermineGapHelper(gapZhistpos,gapZfitslopepos,gapZfitinterceptpos,gapZdivslopepos);
01064   DetermineGapHelper(gapZhist,gapZfitslope,gapZfitintercept,gapZdivslope);
01065 }
01066 //_____________________________________________________________________________
01067 void StSpaceChargeEbyEMaker::DetermineGapHelper(TH2F* gh,
01068       float& fitslope, float& fitintercept, float& divslope) {
01069 
01070   ga1.SetParameters(gh->GetEntries()/(16.*2.*10.),0.,0.1);
01071   ga1.SetParLimits(0,0.001,10.0*gh->GetEntries()); // Loglikelihood only works with positive functions
01072   gh->FitSlicesX(&ga1,1,0,20,"L0Q"); // gapZhist bin contents are integers
01073   const char* hn = gh->GetName();
01074   TH1D* GapsChi2 = (TH1D*) gDirectory->Get(Form("%s_chi2",hn));
01075   TH1D* GapsAmp  = (TH1D*) gDirectory->Get(Form("%s_0",hn));
01076   TH1D* GapsMean = (TH1D*) gDirectory->Get(Form("%s_1",hn));
01077   TH1D* GapsRMS  = (TH1D*) gDirectory->Get(Form("%s_2",hn));
01078 
01079   divslope = GapsMean->GetBinContent(GZN);
01080   egapZdivslope = TMath::Abs(GapsMean->GetBinError(GZN));
01081 
01082   // Fit only from |Z| = 25-150cm
01083   // Zero bins -25 thru 25cm (depends upon GZN,GZL,GZH definitions)
01084   GapsMean->SetBinContent(8,0); GapsMean->SetBinError(8,0);
01085   GapsMean->SetBinContent(9,0); GapsMean->SetBinError(9,0);
01086 
01087   GapsMean->Fit(&ln1,"R0Q");
01088   fitslope = ln1.GetParameter(1);
01089   egapZfitslope = TMath::Abs(ln1.GetParError(1));
01090   fitintercept = ln1.GetParameter(0);
01091   egapZfitintercept = TMath::Abs(ln1.GetParError(0));
01092 
01093   delete GapsChi2;
01094   delete GapsAmp;
01095   delete GapsMean;
01096   delete GapsRMS;
01097 }
01098 //_____________________________________________________________________________
01099 float StSpaceChargeEbyEMaker::EvalCalib(TDirectory* hdir) {
01100 
01101   if (hdir) {
01102     dcehist = (TH2F*) (hdir->Get("DcaEve"));
01103     timehist = (TH1F*) (hdir->Get("EvtTime"));
01104     scehist = (TH1F*) (hdir->Get("SpcChgEvt"));
01105     if (!(dcehist && timehist && scehist)) {
01106       LOG_ERROR << "Problems finding SC histograms" << endm;
01107       return 999.;
01108     }
01109   }
01110 
01111   TH1D* dcaproj = dcehist->ProjectionY();
01112 
01113   // Initial fits to DCA distribution
01114   ga1.SetParameters(1.,0.,1.);
01115   dcaproj->Fit(&ga1,"Q");
01116   float hm = ga1.GetParameter(1);
01117   float hw = TMath::Abs(ga1.GetParameter(2));
01118   float hd = 0.6*hw;
01119   ga1.SetRange(hm-hd,hm+hd);
01120   dcaproj->Fit(&ga1,"RQ");
01121   float gm = ga1.GetParameter(1);
01122   float gw = TMath::Abs(ga1.GetParameter(2));
01123   float gm1 = gm;
01124   float gw1 = gw;
01125 
01126   // Iterate fit to get best answer
01127   ga1.SetRange(gm-0.9*gw,gm+0.9*gw);
01128   dcaproj->Fit(&ga1,"RQ");
01129   gm = ga1.GetParameter(1);
01130   gw = TMath::Abs(ga1.GetParameter(2));
01131   ga1.SetRange(gm-0.8*gw,gm+0.8*gw);
01132   dcaproj->Fit(&ga1,"RQ");
01133   gm = ga1.GetParameter(1);
01134   gw = TMath::Abs(ga1.GetParameter(2));
01135   ga1.SetRange(gm-0.7*gw,gm+0.7*gw);
01136   dcaproj->Fit(&ga1,"RQ");
01137   gm = ga1.GetParameter(1);
01138   gw = TMath::Abs(ga1.GetParameter(2));
01139   float gme = TMath::Abs(ga1.GetParError(1));
01140   float gwe = TMath::Abs(ga1.GetParError(2));
01141 
01142   // Average SC
01143   scehist->Fit("pol0","Q");
01144   TF1* pl0 = scehist->GetFunction("pol0");
01145   float pm = 0, pw = 0;
01146   if (pl0) {
01147     pm = pl0->GetParameter(0);
01148     pw = pl0->GetParError(0);
01149   }
01150   // Other counts
01151   float spc = (float) (scehist->GetEntries());
01152   float dcc = (float) (dcehist->GetEntries());
01153   float evc = (float) (timehist->GetEntries());
01154   timehist->GetXaxis()->SetRange(1,(int) evc);
01155   timehist->Fit("pol0","LQ");
01156   pl0 = timehist->GetFunction("pol0");
01157   float epsec = 0;
01158   if (pl0) epsec = pl0->GetParameter(0); // events per second
01159 
01160   // Quality measures
01161   float wid = TMath::Min(10.,TMath::Log10(gw1*gw1+gw*gw));
01162   float frac = spc/evc; // fraction of events for which SC was found
01163 
01164   float code=0;
01165   if (frac<0.2) code = 1. + frac;               // code = 1.x
01166   else if (wid>0) code = 2. + 0.1*wid;          // code = 2.x
01167 
01168   LOG_INFO << Form("SCeval: %f %f %f %f %f %f %f %f %f %f %f %f %f %f\n",
01169     hm,hw,gm,gw,pm,pw,gm1,gw1,spc,dcc,evc,epsec,gme,gwe) << endm;
01170 
01171   if (code>0) {
01172     LOG_ERROR << "CheckFail: Break of SpaceCharge performance! code = " << code << endm;
01173   }
01174 
01175   return code;
01176 }
01177 //_____________________________________________________________________________
01178 // $Id: StSpaceChargeEbyEMaker.cxx,v 1.41 2012/04/25 19:23:55 genevb Exp $
01179 // $Log: StSpaceChargeEbyEMaker.cxx,v $
01180 // Revision 1.41  2012/04/25 19:23:55  genevb
01181 // Pointing angle near vertex needed for improved PredictSpaceCharge
01182 //
01183 // Revision 1.40  2012/01/14 00:21:04  genevb
01184 // Add code for EMC match, set default to required
01185 //
01186 // Revision 1.39  2011/10/27 23:11:03  genevb
01187 // Account for pointing error in sDCA and predicted SpaceCharge
01188 //
01189 // Revision 1.38  2011/10/26 15:33:49  genevb
01190 // Avoid floating point exception in Loglikelihood fits of gaussian peaks
01191 //
01192 // Revision 1.37  2011/04/18 17:36:52  genevb
01193 // Expanded range for sc hists
01194 //
01195 // Revision 1.36  2011/02/10 18:31:45  genevb
01196 // Restore corrected coincidence rates, add QA histogram of where events/tracks are cut
01197 //
01198 // Revision 1.35  2011/02/09 21:56:50  genevb
01199 // Version which can work in SL10k
01200 //
01201 // Revision 1.34  2011/02/09 21:11:36  genevb
01202 // Parameters need to be available in normal event-by-event mode too
01203 //
01204 // Revision 1.33  2011/02/09 16:24:18  genevb
01205 // Allow for historical operating parameters in Prepass mode
01206 //
01207 // Revision 1.32  2010/11/17 17:23:33  genevb
01208 // Include corrected coincidence rates in ntuple
01209 //
01210 // Revision 1.31  2010/07/09 19:01:54  genevb
01211 // Add TOF matching
01212 //
01213 // Revision 1.30  2010/06/09 20:24:53  genevb
01214 // Modify interface to allow EMC and TOF matching requirements (needs implementation)
01215 //
01216 // Revision 1.29  2010/02/25 21:50:18  genevb
01217 // Pass sector number to StMagUtilities, correct unsigned int usage
01218 //
01219 // Revision 1.28  2010/02/24 22:58:21  genevb
01220 // Even the highest ranked vertex might need to be above a minimum ranking
01221 //
01222 // Revision 1.27  2010/02/23 23:59:41  genevb
01223 // Reduce positional biases in GridLeak measurements
01224 //
01225 // Revision 1.26  2010/01/27 15:11:00  fisyak
01226 // eliminate access to StTpcDbMaker, use directly gStTpcDb
01227 //
01228 // Revision 1.25  2010/01/27 14:42:33  fisyak
01229 // Use  StMagUtilities::Instance instead of asking tpcHitMoverMaker
01230 //
01231 // Revision 1.24  2009/11/16 22:02:19  genevb
01232 // Loosen nDaughters cut, add BEMCmatch cut, PCT hits cut, enable padrow 13 for Run 9+
01233 //
01234 // Revision 1.23  2009/11/10 20:54:13  fisyak
01235 // pams Cleanup
01236 //
01237 // Revision 1.22  2008/07/24 19:17:55  genevb
01238 // SpaceChargeEWRatio must be written to prepass output table
01239 //
01240 // Revision 1.21  2008/07/21 21:17:10  genevb
01241 // No call to EvalCalib() if Prepass runs successfully
01242 //
01243 // Revision 1.20  2008/07/17 23:35:33  jeromel
01244 // Small format change
01245 //
01246 // Revision 1.19  2008/07/16 03:42:15  genevb
01247 // Use CheckFail in failed performance check error message
01248 //
01249 // Revision 1.18  2008/07/15 22:30:38  genevb
01250 // Added evaluation of calibration performance
01251 //
01252 // Revision 1.17  2008/04/30 14:52:15  genevb
01253 // Reduce pileup contributions
01254 //
01255 // Revision 1.16  2008/01/14 19:22:49  genevb
01256 // Fine tuning of parameters, removal of excess text in messages
01257 //
01258 // Revision 1.15  2007/01/25 19:04:04  perev
01259 // GMT fix
01260 //
01261 // Revision 1.14  2007/01/24 21:42:22  perev
01262 // GMT conversion fixed
01263 //
01264 // Revision 1.13  2006/12/16 01:00:58  genevb
01265 // Better handling of zero magnetic field
01266 //
01267 // Revision 1.12  2006/08/15 23:40:59  genevb
01268 // Averaging was done improperly in DontReset mode
01269 //
01270 // Revision 1.11  2006/07/17 20:13:08  genevb
01271 // Disallow SVT points on tracks
01272 //
01273 // Revision 1.10  2006/07/02 23:22:36  genevb
01274 // Allow for SVT/SSD hits on tracks (necessary for ITTF)
01275 //
01276 // Revision 1.9  2006/06/01 17:27:11  genevb
01277 // Bug fix: gapd and gapf backwards; Improvements: gap fit intercepts, hist and fit ranges
01278 //
01279 // Revision 1.8  2006/01/05 19:12:53  genevb
01280 // Added calib mode
01281 //
01282 // Revision 1.7  2005/12/07 19:45:46  perev
01283 // Histos diconnected from directory. EndCrashFix
01284 //
01285 // Revision 1.6  2005/04/21 19:38:20  genevb
01286 // Additional code for studying SpaceCharge
01287 //
01288 // Revision 1.5  2005/02/16 17:18:06  genevb
01289 // Fill StEvent info on SpaceCharge
01290 //
01291 // Revision 1.4  2004/08/13 20:49:12  genevb
01292 // Improve upon keeping method locked on for each event, and timestamp change
01293 //
01294 // Revision 1.3  2004/08/02 01:19:27  genevb
01295 // minor fixes for getting directories correct
01296 //
01297 // Revision 1.2  2004/07/01 01:46:04  genevb
01298 // Slightly larger margin for full vs half field differences
01299 //
01300 // Revision 1.1  2004/06/30 23:16:00  genevb
01301 // Introduction of StSpaceChargeEbyEMaker
01302 //
01303 //