// $Id$ // // File: JEventProcessor_GX_Millepede_Alignment.cc // Created: Mon Apr 20 13:16:44 EDT 2015 // Creator: mstaib (on Linux ifarm1401 2.6.32-431.el6.x86_64 x86_64) // #include "JEventProcessor_GX_Millepede_Alignment.h" using namespace jana; // Routine used to create our JEventProcessor #include #include extern "C"{ void InitPlugin(JApplication *app){ InitJANAPlugin(app); app->AddProcessor(new JEventProcessor_GX_Millepede_Alignment()); } } // "C" //------------------ // JEventProcessor_GX_Millepede_Alignment (Constructor) //------------------ JEventProcessor_GX_Millepede_Alignment::JEventProcessor_GX_Millepede_Alignment() { } //------------------ // ~JEventProcessor_GX_Millepede_Alignment (Destructor) //------------------ JEventProcessor_GX_Millepede_Alignment::~JEventProcessor_GX_Millepede_Alignment() { } //------------------ // init //------------------ jerror_t JEventProcessor_GX_Millepede_Alignment::init(void) { // This is called once at program startup. If you are creating // and filling historgrams in this plugin, you should lock the // ROOT mutex like this: // // japp->RootWriteLock(); // ... fill historgrams or trees ... // japp->RootUnLock(); // milleWriter = new Mille("mille_out.mil"); return NOERROR; } //------------------ // brun //------------------ jerror_t JEventProcessor_GX_Millepede_Alignment::brun(JEventLoop *eventLoop, int runnumber) { // This is called whenever the run number changes map cdc_res_parms; eventLoop->GetCalib("CDC/cdc_resolution_parms", cdc_res_parms); CDC_RES_PAR1 = cdc_res_parms["res_par1"]; CDC_RES_PAR2 = cdc_res_parms["res_par2"]; vector< map > tvals; cdc_drift_table.clear(); if (eventLoop->GetCalib("CDC/cdc_drift_table", tvals)==false){ for(unsigned int i=0; i &row = tvals[i]; cdc_drift_table.push_back(1000.*row["t"]); } } // This is called whenever the run number changes DApplication* dapp=dynamic_cast(eventLoop->GetJApplication()); dgeom = dapp->GetDGeometry(runnumber); //bfield = dapp->GetBfield(); // Get the position of the CDC downstream endplate from DGeometry //double endplate_z,endplate_dz,endplate_rmin,endplate_rmax; //dgeom->GetCDCEndplate(endplate_z,endplate_dz,endplate_rmin,endplate_rmax); //dgeom->GetCDCWires (cdcwires); dgeom->GetCDCCenterZ(cdc_center_z); ///< z-location of center of CDC wires in cm dgeom->GetCDCAxialLength(cdc_axial_length); return NOERROR; } //------------------ // evnt //------------------ jerror_t JEventProcessor_GX_Millepede_Alignment::evnt(JEventLoop *loop, int eventnumber) { // This is called for every event. Use of common resources like writing // to a file or filling a histogram should be mutex protected. Using // loop->Get(...) to get reconstructed objects (and thereby activating the // reconstruction algorithm) should be done outside of any mutex lock // since multiple threads may call this method at the same time. // Here's an example: // // vector mydataclasses; // loop->Get(mydataclasses); // // japp->RootWriteLock(); // ... fill historgrams or trees ... // japp->RootUnLock(); // This should get the straight line tracks in the case that we have 0 field vector trackVector; loop->Get(trackVector); // Now make a few quality cuts and histogram a few things Fill1DHistogram("Millepede", "Tracks", "Number of Tracks", trackVector.size(), "Number of tracks; Number of Tracks; Events", 11, -0.5, 10.5); // Only interested in single track event for now // if (trackVector.size() != 1) return NOERROR; //Only fill these labels once int label[14088]; for (int temp = 0; temp < 14088; temp++){ label[temp]=temp+1; } for (unsigned int i = 0; i < trackVector.size(); i++){ const DTrackTimeBased *thisTrack = trackVector[i]; double FOM = TMath::Prob(thisTrack->chisq, thisTrack->Ndof); Fill1DHistogram("Millepede", "Tracks", "Tracking FOM", FOM, "Tracking FOM; FOM; Entries", 100, 0, 1.0); // Only use tracks above a certain FOM for the study if( FOM < 1.E-10) continue; //Now just fill Mille structure. Loop over the pulls to go hit by hit through the fit // These arrays are holders for some of the data vector pullVector = thisTrack->pulls; Fill1DHistogram("Millepede", "Tracks", "Pulls Size", pullVector.size(), "Number of Pulls;Number; Entries", 57, -0.5, 56.5); // Get the track parameters for this fit float xslope = thisTrack->px() / thisTrack->pz(); float yslope = thisTrack->py() / thisTrack->pz(); float z = thisTrack->z(); float x0 = thisTrack->x() - xslope*z; float y0 = thisTrack->y() - yslope*z; if (xslope != xslope || yslope != yslope || z != z || x0 != x0 || y0 != y0) continue; // Move on to the next track if any NaN int numCDCHits = 0; for ( unsigned int j=0; j < pullVector.size(); j++){ DTrackFitter::pull_t thisPull = pullVector[j]; if (thisPull.cdc_hit == NULL) continue; //Only interested in CDC hits for now double residual = CalculateResidual(thisTrack, thisPull.cdc_hit->wire, thisPull.d); double error = CalculateError(thisPull.tdrift); Fill1DHistogram("Millepede","Tracks","tdrift from Pulls", thisPull.tdrift, "tdrift from Pulls", 100, -100, 1000); Fill2DHistogram("Millepede", "Tracks", "Residual Vs. Drift Distance", thisPull.d, thisPull.resi, "Residual Vs Drift Distance; Drift Distance [cm]; Residual (10 #mum/bin)", 80, 0.0, 0.8, 1000, -0.5, 0.5); Fill2DHistogram("Millepede", "Tracks", "Residual Vs. Drift Time", thisPull.tdrift, thisPull.resi, "Residual Vs Drift Time; Drift Time (2 ns/bin); Residual (10 #mum/bin)", 400, 0.0, 800, 1000, -0.5, 0.5); Fill1DHistogram("Millepede","Tracks","Errors from Calculation", error, "Errors from Calculation", 100, 0.0, 0.05); Fill1DHistogram("Millepede","Tracks","Errors from Pulls", thisPull.err, "Errors from Pulls", 100, 0.0, 0.05); Fill2DHistogram("Millepede", "Tracks", "Residual Vs. Ring Number", thisPull.cdc_hit->wire->ring, residual, "Residual For Each Ring Number; Ring Number; Residual (10 #mum/bin)", 28, 0.5, 28.5, 1000, -0.5, 0.5); // Out to 1cm Fill2DHistogram("Millepede", "Tracks", "Pull Residual Vs. Ring Number", thisPull.cdc_hit->wire->ring, thisPull.resi, "Residual For Each Ring Number; Ring Number; Residual (10 #mum/bin)", 28, 0.5, 28.5, 1000, -0.5, 0.5); // Out to 1cm //if( thisPull.cdc_hit->wire->ring == 28 ) continue; numCDCHits++; int NLC = 4; int NGL = 14088; float localParDer[4] = {}; float globalParDer[14088] = {}; if (thisPull.resi != thisPull.resi || thisPull.err != thisPull.err || residual != residual){ // Check the last thing for NaNs numCDCHits = 0; break; } bool anyNaN = CalculateDerivatives(localParDer, globalParDer, x0, xslope, y0, yslope, thisPull.cdc_hit->wire); milleWriter->mille(NLC, localParDer, NGL, globalParDer, label, residual, thisPull.err); if (anyNaN) { //Notice if we have any NaN values and get out numCDCHits = 0; break; } } Fill1DHistogram("Millepede", "Tracks", "Number of CDC Hits", numCDCHits, "Number of CDC Hits;Number; Entries", 57, -0.5, 56.5); if(numCDCHits > 15) { // If the event is selected, Fill some diagnostic histograms // Need a second loop over the pulls for ( unsigned int j=0; j < pullVector.size(); j++){ DTrackFitter::pull_t thisPull = pullVector[j]; if (thisPull.cdc_hit == NULL) continue; //Only interested in CDC hits for now const DCDCWire * thisWire = thisPull.cdc_hit->wire; // For now, we can just look at the errors as a funtion of ring number Fill2DHistogram("Millepede", "Accepted Tracks", "Residual Vs. Ring Number", thisWire->ring, thisPull.resi, "Residual For Each Ring Number; Ring Number; Residual (10 #mum/bin)", 28, 0.5, 28.5, 1000, -0.5, 0.5); // Out to 1cm } milleWriter->end(); } else milleWriter->kill(); } return NOERROR; } //------------------ // erun //------------------ jerror_t JEventProcessor_GX_Millepede_Alignment::erun(void) { // This is called whenever the run number changes, before it is // changed to give you a chance to clean up before processing // events from the next run number. return NOERROR; } //------------------ // fini //------------------ jerror_t JEventProcessor_GX_Millepede_Alignment::fini(void) { // Called before program exit after event processing is finished. delete milleWriter; return NOERROR; } double JEventProcessor_GX_Millepede_Alignment::CalculateResidual(const DTrackTimeBased * thisTrack, const DCDCWire * thisWire, double measurement){ DVector3 trackPosition = thisTrack->position(); DVector3 trackMomentum = thisTrack->momentum(); DVector3 wirePosition = thisWire->origin; DVector3 wireDirection = thisWire->udir; Float_t a = trackMomentum.Dot(trackMomentum); Float_t b = trackMomentum.Dot(wireDirection); Float_t c = wireDirection.Dot(wireDirection); DVector3 w0 = trackPosition - wirePosition; Float_t d = trackMomentum.Dot(w0); Float_t e = wireDirection.Dot(w0); Float_t sc = ((b*e - c*d)/(a*c-b*b)); DVector3 POCAOnTrack = trackPosition + sc * trackMomentum; DVector3 LOCA = w0 + ((b*e - c*d)/(a*c-b*b))*trackMomentum - ((a*e - b*d)/(a*c-b*b))*wireDirection; Float_t DOCA = LOCA.Mag(); double residual = measurement - DOCA; return residual; } double JEventProcessor_GX_Millepede_Alignment::CalculateError(double tdrift){ double sigma = 0.22516; double cutoffTime = 40.0; if (tdrift>0){ if (tdrift>cdc_drift_table[cdc_drift_table.size()-1]){ // Here the drift time is too large. sigma = 0.22516; } sigma = CDC_RES_PAR1/(tdrift+1.)+CDC_RES_PAR2; // This function is very poorly behaved at low drift times // For times less than cutoffTime assume a linear behavior of our function. if( tdrift < cutoffTime ){ double slope = -1.0 * CDC_RES_PAR1 / (( cutoffTime + 1) * (cutoffTime + 1)); sigma = (CDC_RES_PAR1/(cutoffTime+1.)+CDC_RES_PAR2) + slope * (tdrift - cutoffTime); } } else { // Time is negative, or exactly zero, choose position at wire, with error of t=0 hit double slope = -1.0 * CDC_RES_PAR1 / (( cutoffTime + 1) * (cutoffTime + 1)); sigma = (CDC_RES_PAR1/(cutoffTime+1.)+CDC_RES_PAR2) + slope * (0.0 - cutoffTime); } return sigma; } bool JEventProcessor_GX_Millepede_Alignment::CalculateDerivatives(float *localParDer, float *globalParDer, double x0, double xslope, double y0, double yslope, const DCDCWire * wire){ // Get the wire parameters needed to calculate the derivatives from the wire double xu, xd, yu, yd, zu, zd; double wireLength = wire->L; DVector3 wirePosition = wire->origin; DVector3 wireDirection = wire->udir; DVector3 upstreamPosition = wirePosition - (wireLength/2) * wireDirection; DVector3 downstreamPosition = wirePosition + (wireLength/2) * wireDirection; xu = upstreamPosition.x(); yu = upstreamPosition.y(); zu = upstreamPosition.z(); xd = downstreamPosition.x(); yd = downstreamPosition.y(); zd = downstreamPosition.z(); /* if (wire->straw == 1){ cout << "Downstream Position of Ring " << wire->ring << endl; cout << "x = " << xd << " y = " << yd << " phi = " << TMath::ATan(yd/xd) << endl; } if (wire->straw == 1){ cout << "Printing Key Straw for ring " << wire->ring << endl; cout << "Upstream: (" << xu << "," << yu << "," << zu << ")" << endl; cout << "Downstream: (" << xd << "," << yd << "," << zd << ")" << endl; } */ if (xu != xu || yu != yu || zu != zu || xd != xd || yd != yd || zd != zd) return true; // Return immediately if any NaN int straw_offset[29] = {0,0,42,84,138,192,258,324,404,484,577,670,776,882,1005,1128,1263,1398,1544,1690,1848,2006,2176,2346,2528,2710,2907,3104,3313}; int startIndex = straw_offset[wire->ring] + wire->straw - 1; // Calculate this using wire information // Calculate the derivatives // df/dx0 localParDer[0] = ((yd - yu + yslope*(-zd + zu))*sqrt(pow(x0*yd - x0*yu + xslope*y0*zd - x0*yslope*zd - xslope*yu*zd + xu*(y0 - yd + yslope*zd) - xslope*y0*zu + xslope*yd*zu + x0*yslope*zu - xd*(y0 - yu + yslope*zu),2)/ (pow(yd,2) + pow(xslope,2)*pow(yd,2) + pow(xd,2)*(1 + pow(yslope,2)) + pow(xu,2)*(1 + pow(yslope,2)) - 2*yd*yu - 2*pow(xslope,2)*yd*yu + pow(yu,2) + pow(xslope,2)*pow(yu,2) - 2*yd*yslope*zd + 2*yslope*yu*zd + pow(xslope,2)*pow(zd,2) + pow(yslope,2)*pow(zd,2) - 2*xd*(xu*(1 + pow(yslope,2)) + xslope*(yd*yslope - yslope*yu + zd - zu)) + 2*xslope*xu*(yd*yslope - yslope*yu + zd - zu) + 2*yd*yslope*zu - 2*yslope*yu*zu - 2*pow(xslope,2)*zd*zu - 2*pow(yslope,2)*zd*zu + pow(xslope,2)*pow(zu,2) + pow(yslope,2)*pow(zu,2))))/ (x0*yd - x0*yu + xslope*y0*zd - x0*yslope*zd - xslope*yu*zd + xu*(y0 - yd + yslope*zd) - xslope*y0*zu + xslope*yd*zu + x0*yslope*zu - xd*(y0 - yu + yslope*zu)); if (localParDer[0] != localParDer[0]) return true; // df/dxslope localParDer[1] = (2*(-(((zd - zu)*(-((xd - xu)*(xslope*(x0 - xu) + yslope*(y0 - yu) - zu)) - (x0 - xu)*(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu) + 2*xslope*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu)))/ ((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2))) + (-((x0 - xu)*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2))) + (xd - xu)*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu))/ ((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2)) + ((2*xslope*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - 2*(xd - xu)*(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu))* (zd - zu)*(-((xslope*(x0 - xu) + yslope*(y0 - yu) - zu)*(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu)) + (1 + pow(xslope,2) + pow(yslope,2))*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu)))/ pow((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2),2) - ((2*xslope*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - 2*(xd - xu)*(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu))* (-((pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2))*(xslope*(x0 - xu) + yslope*(y0 - yu) - zu)) + (xslope*(xd - xu) + yslope*(yd - yu) + zd - zu)*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu)))/ pow((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2),2))* (-zu - ((zd - zu)*(-((xslope*(x0 - xu) + yslope*(y0 - yu) - zu)*(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu)) + (1 + pow(xslope,2) + pow(yslope,2))*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu)))/ ((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2)) + (-((pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2))*(xslope*(x0 - xu) + yslope*(y0 - yu) - zu)) + (xslope*(xd - xu) + yslope*(yd - yu) + zd - zu)*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu))/ ((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2))) + 2*(-(((xd - xu)*(-((xd - xu)*(xslope*(x0 - xu) + yslope*(y0 - yu) - zu)) - (x0 - xu)*(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu) + 2*xslope*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu)))/ ((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2))) + (xslope*(-((x0 - xu)*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2))) + (xd - xu)*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu)))/ ((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2)) + ((xd - xu)*(2*xslope*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - 2*(xd - xu)*(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu))* (-((xslope*(x0 - xu) + yslope*(y0 - yu) - zu)*(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu)) + (1 + pow(xslope,2) + pow(yslope,2))*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu)))/ pow((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2),2) - (xslope*(2*xslope*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - 2*(xd - xu)*(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu))* (-((pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2))*(xslope*(x0 - xu) + yslope*(y0 - yu) - zu)) + (xslope*(xd - xu) + yslope*(yd - yu) + zd - zu)*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu)))/ pow((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2),2) + (-((pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2))*(xslope*(x0 - xu) + yslope*(y0 - yu) - zu)) + (xslope*(xd - xu) + yslope*(yd - yu) + zd - zu)*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu))/ ((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2)))* (x0 - xu - ((xd - xu)*(-((xslope*(x0 - xu) + yslope*(y0 - yu) - zu)*(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu)) + (1 + pow(xslope,2) + pow(yslope,2))*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu)))/ ((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2)) + (xslope*(-((pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2))*(xslope*(x0 - xu) + yslope*(y0 - yu) - zu)) + (xslope*(xd - xu) + yslope*(yd - yu) + zd - zu)*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu)))/ ((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2))) + 2*(-(((yd - yu)*(-((xd - xu)*(xslope*(x0 - xu) + yslope*(y0 - yu) - zu)) - (x0 - xu)*(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu) + 2*xslope*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu)))/ ((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2))) + (yslope*(-((x0 - xu)*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2))) + (xd - xu)*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu)))/ ((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2)) + ((yd - yu)*(2*xslope*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - 2*(xd - xu)*(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu))* (-((xslope*(x0 - xu) + yslope*(y0 - yu) - zu)*(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu)) + (1 + pow(xslope,2) + pow(yslope,2))*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu)))/ pow((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2),2) - (yslope*(2*xslope*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - 2*(xd - xu)*(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu))* (-((pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2))*(xslope*(x0 - xu) + yslope*(y0 - yu) - zu)) + (xslope*(xd - xu) + yslope*(yd - yu) + zd - zu)*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu)))/ pow((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2),2))* (y0 - yu - ((yd - yu)*(-((xslope*(x0 - xu) + yslope*(y0 - yu) - zu)*(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu)) + (1 + pow(xslope,2) + pow(yslope,2))*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu)))/ ((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2)) + (yslope*(-((pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2))*(xslope*(x0 - xu) + yslope*(y0 - yu) - zu)) + (xslope*(xd - xu) + yslope*(yd - yu) + zd - zu)*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu)))/ ((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2))))/ (2.*sqrt(pow(-zu - ((zd - zu)*(-((xslope*(x0 - xu) + yslope*(y0 - yu) - zu)*(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu)) + (1 + pow(xslope,2) + pow(yslope,2))*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu)))/ ((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2)) + (-((pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2))*(xslope*(x0 - xu) + yslope*(y0 - yu) - zu)) + (xslope*(xd - xu) + yslope*(yd - yu) + zd - zu)*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu))/ ((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2)),2) + pow(x0 - xu - ((xd - xu)*(-((xslope*(x0 - xu) + yslope*(y0 - yu) - zu)*(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu)) + (1 + pow(xslope,2) + pow(yslope,2))*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu)))/ ((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2)) + (xslope*(-((pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2))*(xslope*(x0 - xu) + yslope*(y0 - yu) - zu)) + (xslope*(xd - xu) + yslope*(yd - yu) + zd - zu)*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu)))/ ((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2)),2) + pow(y0 - yu - ((yd - yu)*(-((xslope*(x0 - xu) + yslope*(y0 - yu) - zu)*(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu)) + (1 + pow(xslope,2) + pow(yslope,2))*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu)))/ ((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2)) + (yslope*(-((pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2))*(xslope*(x0 - xu) + yslope*(y0 - yu) - zu)) + (xslope*(xd - xu) + yslope*(yd - yu) + zd - zu)*((x0 - xu)*(xd - xu) + (y0 - yu)*(yd - yu) - (zd - zu)*zu)))/ ((1 + pow(xslope,2) + pow(yslope,2))*(pow(xd - xu,2) + pow(yd - yu,2) + pow(zd - zu,2)) - pow(xslope*(xd - xu) + yslope*(yd - yu) + zd - zu,2)),2))); if (localParDer[1] != localParDer[1]) return true; // df/dy0 localParDer[2] = ((xd - xu + xslope*(-zd + zu))*sqrt(pow(x0*yd - x0*yu + xslope*y0*zd - x0*yslope*zd - xslope*yu*zd + xu*(y0 - yd + yslope*zd) - xslope*y0*zu + xslope*yd*zu + x0*yslope*zu - xd*(y0 - yu + yslope*zu),2)/ ((1 + pow(xslope,2))*pow(yd,2) + pow(xd,2)*(1 + pow(yslope,2)) + pow(xu,2)*(1 + pow(yslope,2)) - 2*yd*yu - 2*pow(xslope,2)*yd*yu + pow(yu,2) + pow(xslope,2)*pow(yu,2) - 2*yd*yslope*zd + 2*yslope*yu*zd + pow(xslope,2)*pow(zd,2) + pow(yslope,2)*pow(zd,2) - 2*xd*(xu*(1 + pow(yslope,2)) + xslope*(yd*yslope - yslope*yu + zd - zu)) + 2*xslope*xu*(yd*yslope - yslope*yu + zd - zu) - 2*(pow(xslope,2)*zd + yslope*(-yd + yu + yslope*zd))*zu + (pow(xslope,2) + pow(yslope,2))*pow(zu,2))))/ (-(x0*yd) + x0*yu - xslope*y0*zd + x0*yslope*zd + xslope*yu*zd - xu*(y0 - yd + yslope*zd) + xslope*y0*zu - xslope*yd*zu - x0*yslope*zu + xd*(y0 - yu + yslope*zu)); if (localParDer[2] != localParDer[2]) return true; //df/dyslope localParDer[3] = ((xd - xu + xslope*(-zd + zu))*(-(x0*yd) + x0*yu - xslope*y0*zd + x0*yslope*zd + xslope*yu*zd - xu*(y0 - yd + yslope*zd) + xslope*y0*zu - xslope*yd*zu - x0*yslope*zu + xd*(y0 - yu + yslope*zu))* (-(pow(xu,2)*y0*yslope) - x0*xu*yd*yslope + pow(xu,2)*yd*yslope + x0*xu*yslope*yu - x0*xu*zd + pow(xu,2)*zd + y0*yd*zd - y0*yu*zd - yd*yu*zd + pow(yu,2)*zd - y0*yslope*pow(zd,2) + yslope*yu*pow(zd,2) + (x0*xu + (-y0 + yd)*(yd - yu) - yslope*(-2*y0 + yd + yu)*zd)*zu + (-y0 + yd)*yslope*pow(zu,2) + pow(xd,2)*(-(y0*yslope) + yslope*yu + zu) + xd*(xu*(2*y0*yslope - yslope*(yd + yu) - zd - zu) + x0*(yd*yslope - yslope*yu + zd - zu) + xslope*((y0 - yu)*(yd - yu) - zd*zu + pow(zu,2))) - xslope*(x0*(pow(yd - yu,2) + pow(zd - zu,2)) + xu*((y0 - yd)*(yd - yu) + zd*(-zd + zu)))))/ (pow((1 + pow(xslope,2))*pow(yd,2) + pow(xd,2)*(1 + pow(yslope,2)) + pow(xu,2)*(1 + pow(yslope,2)) - 2*yd*yu - 2*pow(xslope,2)*yd*yu + pow(yu,2) + pow(xslope,2)*pow(yu,2) - 2*yd*yslope*zd + 2*yslope*yu*zd + pow(xslope,2)*pow(zd,2) + pow(yslope,2)*pow(zd,2) - 2*xd*(xu*(1 + pow(yslope,2)) + xslope*(yd*yslope - yslope*yu + zd - zu)) + 2*xslope*xu*(yd*yslope - yslope*yu + zd - zu) - 2*(pow(xslope,2)*zd + yslope*(-yd + yu + yslope*zd))*zu + (pow(xslope,2) + pow(yslope,2))*pow(zu,2),2)* sqrt(pow(x0*yd - x0*yu + xslope*y0*zd - x0*yslope*zd - xslope*yu*zd + xu*(y0 - yd + yslope*zd) - xslope*y0*zu + xslope*yd*zu + x0*yslope*zu - xd*(y0 - yu + yslope*zu),2)/ ((1 + pow(xslope,2))*pow(yd,2) + pow(xd,2)*(1 + pow(yslope,2)) + pow(xu,2)*(1 + pow(yslope,2)) - 2*yd*yu - 2*pow(xslope,2)*yd*yu + pow(yu,2) + pow(xslope,2)*pow(yu,2) - 2*yd*yslope*zd + 2*yslope*yu*zd + pow(xslope,2)*pow(zd,2) + pow(yslope,2)*pow(zd,2) - 2*xd*(xu*(1 + pow(yslope,2)) + xslope*(yd*yslope - yslope*yu + zd - zu)) + 2*xslope*xu*(yd*yslope - yslope*yu + zd - zu) - 2*(pow(xslope,2)*zd + yslope*(-yd + yu + yslope*zd))*zu + (pow(xslope,2) + pow(yslope,2))*pow(zu,2)))); if (localParDer[3] != localParDer[3]) return true; // df/dxu globalParDer[startIndex*4] = -(((yd - yu + yslope*(-zd + zu))*(x0*yd - x0*yu + xslope*y0*zd - x0*yslope*zd - xslope*yu*zd + xu*(y0 - yd + yslope*zd) - xslope*y0*zu + xslope*yd*zu + x0*yslope*zu - xd*(y0 - yu + yslope*zu))* (-((1 + pow(xslope,2))*(y0 - yd)*yd) - xslope*xu*y0*yslope + xslope*xu*yd*yslope + pow(xd,2)*(1 + pow(yslope,2)) + y0*yu + pow(xslope,2)*y0*yu - yd*yu - pow(xslope,2)*yd*yu + xslope*xu*zd + y0*yslope*zd - 2*yd*yslope*zd + yslope*yu*zd + pow(xslope,2)*pow(zd,2) + pow(yslope,2)*pow(zd,2) + x0*(-(xd*(1 + pow(yslope,2))) + xu*(1 + pow(yslope,2)) + xslope*(yd*yslope - yslope*yu + zd - zu)) - (pow(xslope,2)*zd + yslope*(y0 - yd + yslope*zd))*zu + xd*(-(xu*(1 + pow(yslope,2))) + xslope*(yslope*(y0 - 2*yd + yu) - 2*zd + zu))))/ (pow((1 + pow(xslope,2))*pow(yd,2) + pow(xd,2)*(1 + pow(yslope,2)) + pow(xu,2)*(1 + pow(yslope,2)) - 2*yd*yu - 2*pow(xslope,2)*yd*yu + pow(yu,2) + pow(xslope,2)*pow(yu,2) - 2*yd*yslope*zd + 2*yslope*yu*zd + pow(xslope,2)*pow(zd,2) + pow(yslope,2)*pow(zd,2) - 2*xd*(xu*(1 + pow(yslope,2)) + xslope*(yd*yslope - yslope*yu + zd - zu)) + 2*xslope*xu*(yd*yslope - yslope*yu + zd - zu) - 2*(pow(xslope,2)*zd + yslope*(-yd + yu + yslope*zd))*zu + (pow(xslope,2) + pow(yslope,2))*pow(zu,2),2)* sqrt(pow(x0*yd - x0*yu + xslope*y0*zd - x0*yslope*zd - xslope*yu*zd + xu*(y0 - yd + yslope*zd) - xslope*y0*zu + xslope*yd*zu + x0*yslope*zu - xd*(y0 - yu + yslope*zu),2)/ ((1 + pow(xslope,2))*pow(yd,2) + pow(xd,2)*(1 + pow(yslope,2)) + pow(xu,2)*(1 + pow(yslope,2)) - 2*yd*yu - 2*pow(xslope,2)*yd*yu + pow(yu,2) + pow(xslope,2)*pow(yu,2) - 2*yd*yslope*zd + 2*yslope*yu*zd + pow(xslope,2)*pow(zd,2) + pow(yslope,2)*pow(zd,2) - 2*xd*(xu*(1 + pow(yslope,2)) + xslope*(yd*yslope - yslope*yu + zd - zu)) + 2*xslope*xu*(yd*yslope - yslope*yu + zd - zu) - 2*(pow(xslope,2)*zd + yslope*(-yd + yu + yslope*zd))*zu + (pow(xslope,2) + pow(yslope,2))*pow(zu,2))))); if (globalParDer[startIndex*4] != globalParDer[startIndex*4]) return true; //df/dxd globalParDer[startIndex*4 +1] = ((yd - yu + yslope*(-zd + zu))*(x0*yd - x0*yu + xslope*y0*zd - x0*yslope*zd - xslope*yu*zd + xu*(y0 - yd + yslope*zd) - xslope*y0*zu + xslope*yd*zu + x0*yslope*zu - xd*(y0 - yu + yslope*zu))* (-pow(xu,2) - (1 + pow(xslope,2))*y0*yd - xslope*xu*y0*yslope - xslope*xu*yd*yslope - pow(xu,2)*pow(yslope,2) + y0*yu + pow(xslope,2)*y0*yu + yd*yu + pow(xslope,2)*yd*yu + 2*xslope*xu*yslope*yu - pow(yu,2) - pow(xslope,2)*pow(yu,2) - xslope*xu*zd + y0*yslope*zd - yslope*yu*zd + x0*(-(xd*(1 + pow(yslope,2))) + xu*(1 + pow(yslope,2)) + xslope*(yd*yslope - yslope*yu + zd - zu)) + (2*xslope*xu - yslope*(y0 + yd - 2*yu) + (pow(xslope,2) + pow(yslope,2))*zd)*zu - (pow(xslope,2) + pow(yslope,2))*pow(zu,2) + xd*(xslope*y0*yslope + xu*(1 + pow(yslope,2)) - xslope*(yslope*yu + zu))))/ (pow((1 + pow(xslope,2))*pow(yd,2) + pow(xd,2)*(1 + pow(yslope,2)) + pow(xu,2)*(1 + pow(yslope,2)) - 2*yd*yu - 2*pow(xslope,2)*yd*yu + pow(yu,2) + pow(xslope,2)*pow(yu,2) - 2*yd*yslope*zd + 2*yslope*yu*zd + pow(xslope,2)*pow(zd,2) + pow(yslope,2)*pow(zd,2) - 2*xd*(xu*(1 + pow(yslope,2)) + xslope*(yd*yslope - yslope*yu + zd - zu)) + 2*xslope*xu*(yd*yslope - yslope*yu + zd - zu) - 2*(pow(xslope,2)*zd + yslope*(-yd + yu + yslope*zd))*zu + (pow(xslope,2) + pow(yslope,2))*pow(zu,2),2)* sqrt(pow(x0*yd - x0*yu + xslope*y0*zd - x0*yslope*zd - xslope*yu*zd + xu*(y0 - yd + yslope*zd) - xslope*y0*zu + xslope*yd*zu + x0*yslope*zu - xd*(y0 - yu + yslope*zu),2)/ ((1 + pow(xslope,2))*pow(yd,2) + pow(xd,2)*(1 + pow(yslope,2)) + pow(xu,2)*(1 + pow(yslope,2)) - 2*yd*yu - 2*pow(xslope,2)*yd*yu + pow(yu,2) + pow(xslope,2)*pow(yu,2) - 2*yd*yslope*zd + 2*yslope*yu*zd + pow(xslope,2)*pow(zd,2) + pow(yslope,2)*pow(zd,2) - 2*xd*(xu*(1 + pow(yslope,2)) + xslope*(yd*yslope - yslope*yu + zd - zu)) + 2*xslope*xu*(yd*yslope - yslope*yu + zd - zu) - 2*(pow(xslope,2)*zd + yslope*(-yd + yu + yslope*zd))*zu + (pow(xslope,2) + pow(yslope,2))*pow(zu,2)))); if (globalParDer[startIndex*4+1] != globalParDer[startIndex*4+1]) return true; //df/dyu globalParDer[startIndex*4 + 2] = -(((xd - xu + xslope*(-zd + zu))*(-(x0*yd) + x0*yu - xslope*y0*zd + x0*yslope*zd + xslope*yu*zd - xu*(y0 - yd + yslope*zd) + xslope*y0*zu - xslope*yd*zu - x0*yslope*zu + xd*(y0 - yu + yslope*zu))* (-((1 + pow(xslope,2))*(y0 - yd)*yd) - xslope*xu*y0*yslope + xslope*xu*yd*yslope + pow(xd,2)*(1 + pow(yslope,2)) + y0*yu + pow(xslope,2)*y0*yu - yd*yu - pow(xslope,2)*yd*yu + xslope*xu*zd + y0*yslope*zd - 2*yd*yslope*zd + yslope*yu*zd + pow(xslope,2)*pow(zd,2) + pow(yslope,2)*pow(zd,2) + x0*(-(xd*(1 + pow(yslope,2))) + xu*(1 + pow(yslope,2)) + xslope*(yd*yslope - yslope*yu + zd - zu)) - (pow(xslope,2)*zd + yslope*(y0 - yd + yslope*zd))*zu + xd*(-(xu*(1 + pow(yslope,2))) + xslope*(yslope*(y0 - 2*yd + yu) - 2*zd + zu))))/ (pow((1 + pow(xslope,2))*pow(yd,2) + pow(xd,2)*(1 + pow(yslope,2)) + pow(xu,2)*(1 + pow(yslope,2)) - 2*yd*yu - 2*pow(xslope,2)*yd*yu + pow(yu,2) + pow(xslope,2)*pow(yu,2) - 2*yd*yslope*zd + 2*yslope*yu*zd + pow(xslope,2)*pow(zd,2) + pow(yslope,2)*pow(zd,2) - 2*xd*(xu*(1 + pow(yslope,2)) + xslope*(yd*yslope - yslope*yu + zd - zu)) + 2*xslope*xu*(yd*yslope - yslope*yu + zd - zu) - 2*(pow(xslope,2)*zd + yslope*(-yd + yu + yslope*zd))*zu + (pow(xslope,2) + pow(yslope,2))*pow(zu,2),2)* sqrt(pow(x0*yd - x0*yu + xslope*y0*zd - x0*yslope*zd - xslope*yu*zd + xu*(y0 - yd + yslope*zd) - xslope*y0*zu + xslope*yd*zu + x0*yslope*zu - xd*(y0 - yu + yslope*zu),2)/ ((1 + pow(xslope,2))*pow(yd,2) + pow(xd,2)*(1 + pow(yslope,2)) + pow(xu,2)*(1 + pow(yslope,2)) - 2*yd*yu - 2*pow(xslope,2)*yd*yu + pow(yu,2) + pow(xslope,2)*pow(yu,2) - 2*yd*yslope*zd + 2*yslope*yu*zd + pow(xslope,2)*pow(zd,2) + pow(yslope,2)*pow(zd,2) - 2*xd*(xu*(1 + pow(yslope,2)) + xslope*(yd*yslope - yslope*yu + zd - zu)) + 2*xslope*xu*(yd*yslope - yslope*yu + zd - zu) - 2*(pow(xslope,2)*zd + yslope*(-yd + yu + yslope*zd))*zu + (pow(xslope,2) + pow(yslope,2))*pow(zu,2))))); if (globalParDer[startIndex*4+2] != globalParDer[startIndex*4+2]) return true; //df/dyd globalParDer[startIndex*4 + 3] = ((xd - xu + xslope*(-zd + zu))*pow(pow(x0*yd - x0*yu + xslope*y0*zd - x0*yslope*zd - xslope*yu*zd + xu*(y0 - yd + yslope*zd) - xslope*y0*zu + xslope*yd*zu + x0*yslope*zu - xd*(y0 - yu + yslope*zu),2)/ ((1 + pow(xslope,2))*pow(yd,2) + pow(xd,2)*(1 + pow(yslope,2)) + pow(xu,2)*(1 + pow(yslope,2)) - 2*yd*yu - 2*pow(xslope,2)*yd*yu + pow(yu,2) + pow(xslope,2)*pow(yu,2) - 2*yd*yslope*zd + 2*yslope*yu*zd + pow(xslope,2)*pow(zd,2) + pow(yslope,2)*pow(zd,2) - 2*xd*(xu*(1 + pow(yslope,2)) + xslope*(yd*yslope - yslope*yu + zd - zu)) + 2*xslope*xu*(yd*yslope - yslope*yu + zd - zu) - 2*(pow(xslope,2)*zd + yslope*(-yd + yu + yslope*zd))*zu + (pow(xslope,2) + pow(yslope,2))*pow(zu,2)),1.5)* (-pow(xu,2) - (1 + pow(xslope,2))*y0*yd - xslope*xu*y0*yslope - xslope*xu*yd*yslope - pow(xu,2)*pow(yslope,2) + y0*yu + pow(xslope,2)*y0*yu + yd*yu + pow(xslope,2)*yd*yu + 2*xslope*xu*yslope*yu - pow(yu,2) - pow(xslope,2)*pow(yu,2) - xslope*xu*zd + y0*yslope*zd - yslope*yu*zd + x0*(-(xd*(1 + pow(yslope,2))) + xu*(1 + pow(yslope,2)) + xslope*(yd*yslope - yslope*yu + zd - zu)) + (2*xslope*xu - yslope*(y0 + yd - 2*yu) + (pow(xslope,2) + pow(yslope,2))*zd)*zu - (pow(xslope,2) + pow(yslope,2))*pow(zu,2) + xd*(xslope*y0*yslope + xu*(1 + pow(yslope,2)) - xslope*(yslope*yu + zu))))/ (pow(x0*yd - x0*yu + xslope*y0*zd - x0*yslope*zd - xslope*yu*zd + xu*(y0 - yd + yslope*zd) - xslope*y0*zu + xslope*yd*zu + x0*yslope*zu - xd*(y0 - yu + yslope*zu),2)*(-(x0*yd) + x0*yu - xslope*y0*zd + x0*yslope*zd + xslope*yu*zd - xu*(y0 - yd + yslope*zd) + xslope*y0*zu - xslope*yd*zu - x0*yslope*zu + xd*(y0 - yu + yslope*zu))); if (globalParDer[startIndex*4+3] != globalParDer[startIndex*4+3]) return true; return false; }