Vertex Reconstruction Update

DataModel/FoMCalculator.cpp

@@ -54,7 +54,8 @@ void FoMCalculator::TimePropertiesLnL(double vtxTime, double& vtxFOM)
   // ================
   for( int idigit=0; idigit<this->fVtxGeo->GetNDigits(); idigit++ ){    
     	int detType = this->fVtxGeo->GetDigitType(idigit); 
-      delta = this->fVtxGeo->GetDelta(idigit) - vtxTime;
+      delta=fVtxGeo->GetExtendedResidual(idigit);                  //delta = this->fVtxGeo->GetDelta(idigit)/* - vtxTime*/;
+      //cout<<"delta check: delta, vtxTime, digitdelta, digittime: "<<delta<<", "<<vtxTime<<", "<<fVtxGeo->GetDelta(idigit)<<", "<<fVtxGeo->GetDigitT(idigit) << endl;
       sigma = this->fVtxGeo->GetDeltaSigma(idigit);
       type = this->fVtxGeo->GetDigitType(idigit);
       if (type == RecoDigit::PMT8inch){  //PMT8Inch
@@ -67,6 +68,7 @@ void FoMCalculator::TimePropertiesLnL(double vtxTime, double& vtxFOM)
       }
       A  = 1.0 / ( 2.0*sigma*sqrt(0.5*TMath::Pi()) ); //normalisation constant
       Preal = A*exp(-(delta*delta)/(2.0*sigma*sigma));
+      //cout<<"Pcheck: " << Preal<<endl;
       P = (1.0-Pnoise)*Preal + Pnoise; 
       chi2 += -2.0*log(P);
       ndof += 1.0; 
@@ -80,6 +82,7 @@ void FoMCalculator::TimePropertiesLnL(double vtxTime, double& vtxFOM)
   // return figure of merit
   // ======================
   vtxFOM = fom;
+  //cout<<"timefom: "<< vtxFOM<<endl;
   return;
 }
 
@@ -103,7 +106,7 @@ void FoMCalculator::ConePropertiesFoM(double coneEdge, double& coneFOM)
   double fom = -9999.;
 
   for( int idigit=0; idigit<this->fVtxGeo->GetNDigits(); idigit++ ){ 	
-    if( this->fVtxGeo->IsFiltered(idigit) && this->fVtxGeo->GetDigitType(idigit) == RecoDigit::PMT8inch){
+    if( this->fVtxGeo->IsFiltered(idigit) && this->fVtxGeo->GetDigitType(idigit) == RecoDigit::PMT8inch) {
       deltaAngle = this->fVtxGeo->GetAngle(idigit) - coneEdge;
       digitCharge = this->fVtxGeo->GetDigitQ(idigit);
 
@@ -137,83 +140,6 @@ void FoMCalculator::ConePropertiesFoM(double coneEdge, double& coneFOM)
   return;
 }
 
-void FoMCalculator::ConePropertiesLnL(double vtxX, double vtxY, double vtxZ, double dirX, double dirY, double dirZ, double coneEdge, double& chi2, TH1D angularDist, double& phimax, double& phimin) {
-    double coneEdgeLow = 21.0;  // cone edge (low side)      
-    double coneEdgeHigh = 3.0;  // cone edge (high side)   [muons: 3.0, electrons: 7.0]
-    double deltaAngle = 0.0;
-    double digitCharge = 0.0;
-    double digitPE = 0.0;
-    double coneCharge = 0.0;
-    double allCharge = 0.0;
-    double outerCone = -99.9;
-    double coef = angularDist.Integral(); //1000;
-    chi2 = 0;
-    cout << "ConePropertiesLnL Position: (" << vtxX << ", " << vtxY << ", " << vtxZ << ")" << endl;
-    cout << "And Direction: (" << dirX << ", " << dirY << ", " << dirZ << ")" << endl;
-
-    double digitX, digitY, digitZ;
-    double dx, dy, dz, ds;
-    double px, py, pz;
-    double cosphi, phi, phideg;
-    phimax = 0;
-    phimin = 10;
-    double allPE = 0;
-    int refbin;
-    double weight;
-    double P;
-
-    for (int idigit = 0; idigit < this->fVtxGeo->GetNDigits(); idigit++) {
-        if (this->fVtxGeo->IsFiltered(idigit) && this->fVtxGeo->GetDigitType(idigit) == RecoDigit::PMT8inch) {
-            digitCharge = this->fVtxGeo->GetDigitQ(idigit);
-            allCharge += digitCharge;
-        }
-    }
-
-    for (int idigit = 0; idigit < this->fVtxGeo->GetNDigits(); idigit++) {
-        if (this->fVtxGeo->IsFiltered(idigit) && this->fVtxGeo->GetDigitType(idigit) == RecoDigit::PMT8inch) {
-            deltaAngle = this->fVtxGeo->GetAngle(idigit) - coneEdge;
-            digitCharge = this->fVtxGeo->GetDigitQ(idigit);
-            //digitPE = this->fVtxGeo->GetDigitPE(idigit);
-            digitX = fVtxGeo->GetDigitX(idigit);
-            digitY = fVtxGeo->GetDigitY(idigit);
-            digitZ = fVtxGeo->GetDigitZ(idigit);
-            dx = digitX - vtxX;
-            dy = digitY - vtxY;
-            dz = digitZ - vtxZ;
-            std::cout << "dx, dy, dz: " << dx << ", " << dy << ", " << dz << endl;
-            ds = pow(dx * dx + dy * dy + dz * dz, 0.5);
-            std::cout << "ds: " << ds << endl;
-            px = dx / ds;
-            py = dy / ds;
-            pz = dz / ds;
-            std::cout << "px, py, pz: " << px << ", " << py << ", " << pz << endl;
-            std::cout << "dirX, dirY, DirZ: " << dirX << ", " << dirY << ", " << dirZ << endl;
-
-            cosphi = 1.0;
-            phi = 0.0;
-            //cout << "angle direction: " << dx << " " << dy << " " << dz << " = " << ds << endl;
-            cosphi = px * dirX + py * dirY + pz * dirZ;
-            //cout << "cosphi: " << cosphi << endl;
-            phi = acos(cosphi);
-
-            if (phi > phimax) phimax = phi;
-            if (phi < phimin) phimin = phi;
-
-            phideg = phi / (TMath::Pi() / 180);
-            std::cout << "phi, phideg: " << phi << ", " << phideg << endl;
-            std::cout << "vs. Zenith: " << fVtxGeo->GetZenith(idigit) << endl;
-            refbin = angularDist.FindBin(phideg);
-            weight = angularDist.GetBinContent(refbin)/coef;
-            P = digitCharge / allCharge;
-            //cout << "conefomlnl P: " << P << ", weight: " << weight << endl;
-            chi2 += pow(P - weight, 2)/weight;
-
-            //outerCone = -outhits/inhits;
-        }
-    }
-    //chi2 = (100 - chi2) * exp(-pow(pow(0.7330382, 2) - pow(phimax - phimin, 2), 2) / pow(0.7330382, 2));
-}
-
 
 
 
@@ -369,7 +295,7 @@ void FoMCalculator::PointVertexChi2(double vtxX, double vtxY, double vtxZ,
 
   // calculate residuals
   // ===================
-  this->fVtxGeo->CalcPointResiduals(vtxX, vtxY, vtxZ, 0.0, 
+  this->fVtxGeo->CalcPointResiduals(vtxX, vtxY, vtxZ, vtxTime,
                                  dirX, dirY, dirZ); //calculate expected vertex time for each digit
   // calculate figure of merit
   // =========================
@@ -402,7 +328,7 @@ void FoMCalculator::ExtendedVertexChi2(double vtxX, double vtxY, double vtxZ, do
 
   // calculate residuals
   // ===================
-  this->fVtxGeo->CalcExtendedResiduals(vtxX,vtxY,vtxZ,0.0,dirX,dirY,dirZ);
+  this->fVtxGeo->CalcExtendedResiduals(vtxX,vtxY,vtxZ, vtxTime,dirX,dirY,dirZ);
 
   // calculate figure of merit
   // =========================
@@ -424,39 +350,6 @@ void FoMCalculator::ExtendedVertexChi2(double vtxX, double vtxY, double vtxZ, do
   return;
 }
 
-void FoMCalculator::ExtendedVertexChi2(double vtxX, double vtxY, double vtxZ, double dirX, double dirY, double dirZ, double coneAngle, double vtxTime, double& fom, TH1D pdf)
-{
-	// figure of merit
-	// ===============
-	double vtxFOM = -9999.;
-	double timeFOM = -9999.;
-	double coneFOM = -9999.;
-    double phimax, phimin;
-
-	// calculate residuals
-	// ===================
-	this->fVtxGeo->CalcExtendedResiduals(vtxX, vtxY, vtxZ, 0.0, dirX, dirY, dirZ);
-
-	// calculate figure of merit
-	// =========================
-
-    this->ConePropertiesLnL(vtxX, vtxY, vtxZ, dirX, dirY, dirZ, coneAngle, coneFOM, pdf, phimax, phimin);
-	this->TimePropertiesLnL(vtxTime, timeFOM);
-
-	double fTimeFitWeight = this->fTimeFitWeight;
-	double fConeFitWeight = this->fConeFitWeight;
-	vtxFOM = (fTimeFitWeight*timeFOM + fConeFitWeight * coneFOM) / (fTimeFitWeight + fConeFitWeight);
-
-	// calculate overall figure of merit
-	// =================================
-	fom = vtxFOM;
-
-	// truncate
-	if (fom < -9999.) fom = -9999.;
-
-	return;
-}
-
 
 
 //KEPT FOR HISTORY, BUT FITTER IS CURRENTLY NOT WORKING

DataModel/FoMCalculator.h

@@ -50,7 +50,6 @@ class FoMCalculator {
   double FindSimpleTimeProperties(double myConeEdge);
   void TimePropertiesLnL(double vtxTime, double& vtxFom);
   void ConePropertiesFoM(double coneEdge, double& chi2);
-  void ConePropertiesLnL(double vtxX, double vtxY, double VtxZ, double dirX, double dirY, double dirZ, double coneEdge, double& chi2, TH1D angularDist, double& phimax, double& phimin);
   void PointPositionChi2(double vtxX, double vtxY, double vtxZ, double vtxTime, double& fom);
   void PointDirectionChi2(double vtxX, double vtxY, double vtxZ, double dirX, double dirY, double dirZ, double coneAngle, double& fom);
   void PointVertexChi2(double vtxX, double vtxY, double vtxZ,
@@ -59,9 +58,6 @@ class FoMCalculator {
   void ExtendedVertexChi2(double vtxX, double vtxY, double vtxZ, 
 	                                    double dirX, double dirY, double dirZ, 
 	                                    double coneAngle, double vtxTime, double& fom);
-  void ExtendedVertexChi2(double vtxX, double vtxY, double vtxZ,
-	  double dirX, double dirY, double dirZ,
-	  double coneAngle, double vtxTime, double& fom, TH1D pdf);
 //  void ConePropertiesLnL(double coneParam0, double coneParam1, double coneParam2, double& coneAngle, double& coneFOM);
 //  void CorrectedVertexChi2(double vtxX, double vtxY, double vtxZ, 
 //	                                    double dirX, double dirY, double dirZ, 

DataModel/Hit.h

@@ -15,15 +15,18 @@ class Hit : public SerialisableObject{
 	public:
 	Hit() : TubeId(0), Time(0), Charge(0){serialise=true;}
 	Hit(int thetubeid, double thetime, double thecharge) : TubeId(thetubeid), Time(thetime), Charge(thecharge){serialise=true;}
+	Hit(int hitid, int thetubeid, double thetime, double thecharge) : HitID(hitid),TubeId(thetubeid),Time(thetime),Charge(thecharge){serialise=true;}
 	virtual ~Hit(){};
 
 	inline int GetTubeId() const {return TubeId;}
 	inline double GetTime() const {return Time;}
 	inline double GetCharge() const {return Charge;}
+	inline int GetHitID() const {return HitID;}
 
 	inline void SetTubeId(int tubeid){TubeId=tubeid;}
 	inline void SetTime(double tc){Time=tc;}
 	inline void SetCharge(double chg){Charge=chg;}
+	inline void SetHitID(int hitid){HitID=hitid;}
 
 	bool Print() {
 	  std::cout<<"TubeId : "<<TubeId<<endl;
@@ -33,6 +36,7 @@ class Hit : public SerialisableObject{
 	}
 
 	protected:
+	int HitID;
 	int TubeId;
 	double Time;
 	double Charge;
@@ -59,6 +63,7 @@ class MCHit : public Hit {
 	public:
 	MCHit() : Hit(), Parents(std::vector<int>{}) {serialise=true;}
 	MCHit(int tubeid, double thetime, double thecharge, std::vector<int> theparents) : Hit(tubeid, thetime, thecharge), Parents(theparents) {serialise=true;}
+	MCHit(int hitid, int tubeid, double thetime, double thecharge, std::vector<int> theparents) : Hit(hitid, tubeid, thetime, thecharge), Parents(theparents) { serialise = true; }
 	virtual ~MCHit(){};
 
 	const std::vector<int>* GetParents() const { return &Parents; }

DataModel/MinuitOptimizer.cpp

@@ -904,155 +904,6 @@ void MinuitOptimizer::FitExtendedVertexWithMinuit() {
   return;
 }
 
-void MinuitOptimizer::FitExtendedVertexWithMinuit(TH1D pdf) {
-    // seed vertex
-    // ===========
-    bool foundSeed = (fSeedVtx->FoundVertex()
-        && fSeedVtx->FoundDirection());
-
-
-    double seedTime = fSeedVtx->GetTime();
-    double seedX = fSeedVtx->GetPosition().X();
-    double seedY = fSeedVtx->GetPosition().Y();
-    double seedZ = fSeedVtx->GetPosition().Z();
-
-    double seedDirX = fSeedVtx->GetDirection().X();
-    double seedDirY = fSeedVtx->GetDirection().Y();
-    double seedDirZ = fSeedVtx->GetDirection().Z();
-
-    double seedTheta = acos(seedDirZ);
-    double seedPhi = 0.0;
-
-    //modified by JW
-    if (seedDirX > 0.0) {
-        seedPhi = atan(seedDirY / seedDirX);
-    }
-    if (seedDirX < 0.0) {
-        seedPhi = atan(seedDirY / seedDirX);
-        if (seedDirY > 0.0) seedPhi += TMath::Pi();
-        if (seedDirY <= 0.0) seedPhi -= TMath::Pi();
-    }
-    if (seedDirX == 0.0) {
-        if (seedDirY > 0.0) seedPhi = 0.5 * TMath::Pi();
-        else if (seedDirY < 0.0) seedPhi = -0.5 * TMath::Pi();
-        else seedPhi = 0.0;
-    }
-    // current status
-    // ==============
-    int status = fSeedVtx->GetStatus();
-
-    // reset counter
-    // =============
-    extended_vertex_reset_itr();
-
-    // abort if necessary
-    // ==================
-    if (foundSeed == 0) {
-        if (fPrintLevel >= 0) {
-            std::cout << "   <warning> extended vertex fit failed to find input vertex " << std::endl;
-        }
-        status |= RecoVertex::kFailExtendedVertex;
-        fFittedVtx->SetStatus(status);
-        return;
-    }
-
-    // run Minuit
-    // ==========  
-    // six-parameter fit to vertex position, time and direction
-
-    int err = 0;
-    int flag = 0;
-
-    double fitXpos = 0.0;
-    double fitYpos = 0.0;
-    double fitZpos = 0.0;
-    double fitTheta = 0.0;
-    double fitPhi = 0.0;
-    double fitTime = 0.0;
-
-    double fitTimeErr = 0.0;
-    double fitXposErr = 0.0;
-    double fitYposErr = 0.0;
-    double fitZposErr = 0.0;
-    double fitThetaErr = 0.0;
-    double fitPhiErr = 0.0;
-
-    double* arglist = new double[10];
-    arglist[0] = 2;  // 1: standard minimization
-    // 2: try to improve minimum
-
-// re-initialize everything...
-    fMinuitExtendedVertex->mncler();
-    fMinuitExtendedVertex->SetFCN(extended_vertex_chi2);
-    fMinuitExtendedVertex->mnset();
-    fMinuitExtendedVertex->mnexcm("SET STR", arglist, 1, err);
-    fMinuitExtendedVertex->mnparm(0, "x", seedX, 1.0, fXmin, fXmax, err);
-    fMinuitExtendedVertex->mnparm(1, "y", seedY, 1.0, fYmin, fYmax, err);
-    fMinuitExtendedVertex->mnparm(2, "z", seedZ, 5.0, fZmin, fZmax, err);
-    fMinuitExtendedVertex->mnparm(3, "theta", seedTheta, 0.125 * TMath::Pi(), -1.0 * TMath::Pi(), 2.0 * TMath::Pi(), err);
-    fMinuitExtendedVertex->mnparm(4, "phi", seedPhi, 0.125 * TMath::Pi(), -2.0 * TMath::Pi(), 2.0 * TMath::Pi(), err);
-    fMinuitExtendedVertex->mnparm(5, "vtxTime", seedTime, 1.0, fTmin, fTmax, err); //....TX
-
-    flag = fMinuitExtendedVertex->Migrad();
-
-    fMinuitExtendedVertex->GetParameter(0, fitXpos, fitXposErr);
-    fMinuitExtendedVertex->GetParameter(1, fitYpos, fitYposErr);
-    fMinuitExtendedVertex->GetParameter(2, fitZpos, fitZposErr);
-    fMinuitExtendedVertex->GetParameter(3, fitTheta, fitThetaErr);
-    fMinuitExtendedVertex->GetParameter(4, fitPhi, fitPhiErr);
-    fMinuitExtendedVertex->GetParameter(5, fitTime, fitTimeErr);
-
-    delete[] arglist;
-
-    //correct angles, JW
-    if (fitTheta < 0.0) fitTheta = -1.0 * fitTheta;
-    if (fitTheta > TMath::Pi()) fitTheta = 2.0 * TMath::Pi() - fitTheta;
-    if (fitPhi < -1.0 * TMath::Pi()) fitPhi += 2.0 * TMath::Pi();
-    if (fitPhi > TMath::Pi()) fitPhi -= 2.0 * TMath::Pi();
-
-    // sort results
-    // ============
-    fVtxX = fitXpos;
-    fVtxY = fitYpos;
-    fVtxZ = fitZpos;
-    fVtxTime = fitTime;
-    fDirX = sin(fitTheta) * cos(fitPhi);
-    fDirY = sin(fitTheta) * sin(fitPhi);
-    fDirZ = cos(fitTheta);
-
-    fVtxFOM = -9999.0;
-
-    fPass = 0;               // flag = 0: normal termination
-    if (flag == 0) fPass = 1; // anything else: abnormal termination 
-
-    fItr = extended_vertex_iterations();
-
-    // fit complete; calculate fit results
-    // ================
-    fgFoMCalculator->ExtendedVertexChi2(fVtxX, fVtxY, fVtxZ,
-        fDirX, fDirY, fDirZ,
-        fConeAngle, fVtxTime, fVtxFOM, pdf);
-
-    // set vertex and direction
-    // ========================
-    fFittedVtx->SetVertex(fVtxX, fVtxY, fVtxZ, fVtxTime);
-    fFittedVtx->SetDirection(fDirX, fDirY, fDirZ);
-    fFittedVtx->SetConeAngle(fConeAngle);
-    fFittedVtx->SetFOM(fVtxFOM, fItr, fPass);
-
-    // set status
-    // ==========
-    bool inside_det = ANNIEGeometry::Instance()->InsideDetector(fVtxX, fVtxY, fVtxZ);
-    if (!fPass || !inside_det) status |= RecoVertex::kFailExtendedVertex;
-    fFittedVtx->SetStatus(status);
-    if (fPrintLevel >= 0) {
-        if (!fPass) std::cout << "   <warning> extended vertex fit failed to converge! Error code: " << flag << std::endl;
-    }
-    // return vertex
-    // =============  
-    return;
-}
-
 
 //KEPT FOR HISTORY, BUT FITTER IS CURRENTLY NOT WORKING
 //THESE SHOULD BE MOVED TO BEFORE THE CONSTRUCTOR

DataModel/MinuitOptimizer.h

@@ -87,7 +87,6 @@ class MinuitOptimizer {
   void FitPointDirectionWithMinuit();
   void FitPointVertexWithMinuit();
   void FitExtendedVertexWithMinuit();
-  void FitExtendedVertexWithMinuit(TH1D pdf);
 
   double GetTime() {return fVtxTime;}
   double GetFOM() {return fVtxFOM;}