// $Id$ // // File: DParticle_factory_THROWN.cc // Created: Sat Oct 4 22:04:56 EDT 2008 // Creator: davidl (on Darwin Amelia.local 8.11.1 i386) // #include using namespace std; #include #include #include #include "DParticle_factory_THROWN.h" #include #include #include #include //------------------ // Constructor //------------------ DParticle_factory_THROWN::DParticle_factory_THROWN() { DEBUG_LEVEL = 0; // Define target center target = new DCoordinateSystem(); target->origin.SetXYZ(0.0, 0.0, 65.0); target->sdir.SetXYZ(1.0, 0.0, 0.0); target->tdir.SetXYZ(0.0, 1.0, 0.0); target->udir.SetXYZ(0.0, 0.0, 1.0); target->L = 30.0; } //------------------ // brun //------------------ jerror_t DParticle_factory_THROWN::brun(jana::JEventLoop *eventLoop, int runnumber) { // Get resolutions of chambers. These should be set somewhere in the TRACKING package SIGMA_CDC = 0.0; SIGMA_FDC_ANODE = 0.0; SIGMA_FDC_CATHODE = 0.0; JParameter *parm = gPARMS->GetParameter("TRKFIT:SIGMA_CDC"); if(parm)SIGMA_CDC = parm->d(); parm = gPARMS->GetParameter("TRKFIT:SIGMA_FDC_ANODE"); if(parm)SIGMA_FDC_ANODE = parm->d(); parm = gPARMS->GetParameter("TRKFIT:SIGMA_FDC_CATHODE"); if(parm)SIGMA_FDC_CATHODE = parm->d(); _DBG_<<" SIGMA_CDC = "<AddAssociatedObject(cdchits[i]); for(unsigned int i=0; iAddAssociatedObject(fdchits[i]); // Calculate chisq and Ndof vector wires; vector shifts; vector errs; GetWiresShiftsErrs(kTimeBased, rt, wires, shifts, errs); vector chisqv; double chisq; ChiSq(rt, wires, shifts, errs, chisqv, &chisq, &particle->Ndof); particle->chisq = chisq; _data.push_back(particle); } return NOERROR; } //============================================================== // NOTE: The following routines were copied from DParticle_factory //============================================================== //------------------ // AddCDCTrackHits //------------------ void DParticle_factory_THROWN::AddCDCTrackHits(DReferenceTrajectory *rt, vector &cdctrackhits) { /// Determine the probability that for each CDC hit that it came from the particle with the given trajectory. /// /// This will calculate a probability for each CDC hit that /// it came from the particle represented by the given /// DReference trajectory. The probability is based on /// the residual between the distance of closest approach /// of the trajectory to the wire and the drift time. // Calculate beta of particle assuming its a pion for now. If the // particles is really a proton or an electron, the residual // calculated below will only be off by a little. double TOF_MASS = 0.13957018; double beta = 1.0/sqrt(1.0+TOF_MASS*TOF_MASS/rt->swim_steps[0].mom.Mag2()); // The error on the residual. This should include the // error from measurement,particle parameters, and multiple // scattering. //double sigma = sqrt(pow(SIGMA_CDC,2.0) + pow(0.4000,2.0)); double sigma = 0.8/sqrt(12.0); // Minimum probability of hit belonging to wire and still be accepted double MIN_HIT_PROB = 0.2; for(unsigned int j=0; jDistToRT(hit->wire, &s); // Distance using drift time // NOTE: Right now we assume pions for the TOF // and a constant drift velocity of 55um/ns double tof = s/(beta*3E10*1E-9); double dist = (hit->tdrift - tof)*55E-4; // Residual //double resi = dist - doca; double resi = doca - 0.4; double chisq = pow(resi/sigma, 2.0); // Use chi-sq probaility function with Ndof=1 to calculate probability double probability = TMath::Prob(chisq/4.0, 1); if(probability>=MIN_HIT_PROB)cdchits.push_back(hit); if(debug_level>10)_DBG_<<"s="< &wires, vector &shifts, vector &errs) { swim_step_t &start_step = rt->swim_steps[0]; DVector3 &pos = start_step.origin; #if 1 // --- Target --- // Define the target as though it were a wire so it is included // in the chi-sq target->L = 3.0; //target->L = 0.0; target->origin.SetXYZ(0.0, 0.0, pos.Z()); wires.push_back(target); errs.push_back(fit_type==kWireBased ? 0.1:0.1); //errs.push_back(fit_type==kWireBased ? 0.0001:0.0001); if(fit_type!=kWireBased)shifts.push_back(DVector3(0.0, 0.0, 0.0)); #endif // --- CDC --- for(unsigned int i=0; iwire; wires.push_back(wire); // Fill in shifts and errs vectors based on whether we're doing // hit-based or time-based tracking if(fit_type==kWireBased){ // If we're doing hit-based tracking then only the wire positions // are used and the drift time info is ignored. Thus, we don't // have to calculate the shift vectors //errs.push_back(0.261694); // empirically from single pi+ events errs.push_back(0.8/sqrt(12.0)); // variance for evenly illuminated straw }else{ // Find the DOCA point for the RT and use the momentum direction // there and the wire direction to determine the "shift". // Note that whether the shift is to the left or to the right // is not decided here. That ambiguity is left to be resolved later // by applying a minus sign (or not) to the shift. DVector3 pos_doca, mom_doca; double s; rt->DistToRT(wire, &s); rt->GetLastDOCAPoint(pos_doca, mom_doca); DVector3 shift = wire->udir.Cross(mom_doca); // The magnitude of the shift is based on the drift time. The // value of the dist member of the DCDCTrackHit object does not // subtract out the TOF. This can add 50-100 microns to the // resolution in the CDC. Here, we actually can calculate the TOF // (for a given mass hypothesis). double mass = 0.13957; double beta = 1.0/sqrt(1.0 + pow(mass/mom_doca.Mag(), 2.0))*2.998E10; double tof = s/beta/1.0E-9; // in ns double dist = hit->dist*((hit->tdrift-tof)/hit->tdrift); shift.SetMag(dist); // If we're doing time-based tracking then there is a sign ambiguity // to each of the "shifts". It is not realistic to check every possible // solution so we have to make a guess as to what the sign is for each // based on the current reference trajectory. Presumably, if, we're // doing time-based track fitting then we have already gone through a pass // of hit-based track fitting so the initial reference trajectory should be // a pretty good indicator as to what side of the wire the track went // on. We use this in the assignments for now. In the future, we should // try multiple cases for hits close to the wire where the ambiguity // is not clearly resolved. // shift needs to be in the direction pointing from the avalanche // position on the wire to the DOCA point on the rt. We find this // by first finding this vector for the rt and then seeing if // the "shift" vector is generally pointing in the same direction // as it or not by checking if thier dot product is positive or // negative. double u = rt->GetLastDistAlongWire(); DVector3 pos_wire = wire->origin + u*wire->udir; DVector3 pos_diff = pos_doca-pos_wire; if(shift.Dot(pos_diff)<0.0)shift = -shift; shifts.push_back(shift); errs.push_back(SIGMA_CDC); } } // --- FDC --- for(unsigned int i=0; iwire; wires.push_back(wire); // Fill in shifts and errs vectors based on whether we're doing // hit-based or time-based tracking if(fit_type==kWireBased){ // If we're doing hit-based tracking then only the wire positions // are used and the drift time info is ignored. Thus, we don't // have to calculate the shift vectors //errs.push_back(0.261694); // empirically from single pi+ events errs.push_back(0.5/sqrt(12.0)); // variance for evenly illuminated cell }else{ // Find the DOCA point for the RT and use the momentum direction // there and the wire direction to determine the "shift". // Note that whether the shift is to the left or to the right // is not decided here. That ambiguity is left to be resolved later // by applying a minus sign (or not) to the shift. DVector3 pos_doca, mom_doca; double s; rt->DistToRT(wire, &s); rt->GetLastDOCAPoint(pos_doca, mom_doca); DVector3 shift = wire->udir.Cross(mom_doca); // The magnitude of the shift is based on the drift time. The // value of the dist member of the DCDCTrackHit object does not // subtract out the TOF. This can add 50-100 microns to the // resolution in the CDC. Here, we actually can calculate the TOF // (for a given mass hypothesis). double mass = 0.13957; double beta = 1.0/sqrt(1.0 + pow(mass/mom_doca.Mag(), 2.0))*2.998E10; double tof = s/beta/1.0E-9; // in ns double dist = hit->dist*((hit->time-tof)/hit->time); shift.SetMag(dist); // If we're doing time-based tracking then there is a sign ambiguity // to each of the "shifts". It is not realistic to check every possible // solution so we have to make a guess as to what the sign is for each // based on the current reference trajectory. Presumably, if, we're // doing time-based track fitting then we have already gone through a pass // of hit-based track fitting so the initial reference trajectory should be // a pretty good indicator as to what side of the wire the track went // on. We use this in the assignments for now. In the future, we should // try multiple cases for hits close to the wire where the ambiguity // is not clearly resolved. // shift needs to be in the direction pointing from the avalanche // position on the wire to the DOCA point on the rt. We find this // by first finding this vector for the rt and then seeing if // the "shift" vector is generally pointing in the same direction // as it or not by checking if thier dot product is positive or // negative. double u = rt->GetLastDistAlongWire(); DVector3 pos_wire = wire->origin + u*wire->udir; DVector3 pos_diff = pos_doca-pos_wire; if(shift.Dot(pos_diff)<0.0)shift = -shift; shifts.push_back(shift); errs.push_back(SIGMA_FDC_ANODE); } } }