C++ TVector3::SetMagThetaPhi方法代码示例

TGraph GetCurve(int Points,const double & hi_ex_set)
{
  TGraph curve;
  if(!gPrimaryReaction.IsSet()){
    std::cout<<"Reaction Masses have not been set"<<std::endl;
    exit(EXIT_FAILURE);
  }
  if(!gPrimaryReaction.BeamEnergy()){
    std::cout<<"Beam Energy has not been set"<<std::endl;
    exit(EXIT_FAILURE);
  }

  sim::RN_SimEvent evt1(gPrimaryReaction.BeamEnergy(),gPrimaryReaction.M_Beam(),gPrimaryReaction.M_Target(),gPrimaryReaction.M_Recoil(),gPrimaryReaction.M_Fragment()); 
  // Fill the points of the kinematic curve
  int p=0;
  while(p<Points){
    double theta_deg = 180.0*p/Points;
    double phi=2.*M_PI*global::myRnd.Rndm();

    TVector3 nv; nv.SetMagThetaPhi(1.,theta_deg*M_PI/180.0,phi);
    if(!evt1.radiate_in_CM(nv,hi_ex_set))
      continue;
    else 
      curve.SetPoint(p, evt1.getLVrad().Theta()*180/3.14,(double)(evt1.getLVrad().E()-evt1.getLVrad().M()));      
    p++;
  }
  // end for(p)    

  return curve;
}

Double_t getShiftChi2(const Double_t* thetaPhi) {
   Double_t chi2(0), c(0);

   TVector3 norm; norm.SetMagThetaPhi(1.0, thetaPhi[0], thetaPhi[1]);
   for (Int_t ch=0; ch<NSnConstants::kNchans; ++ch) {
      for (Int_t xc=0; xc<ch; ++xc) {
         Double_t dtcor=0;
         for (Int_t i=(ch-xc); i>0; --i) {
            dtcor += dtCorrs[ch-i];
         }
         const TVector3& posCh = getStnPos(ch);
         const TVector3& posXc = getStnPos(xc);
         const Double_t  disCh = -(posCh.Dot(norm));
         const Double_t  disXc = -(posXc.Dot(norm));
         // !!! check sign of delta(distance) and dtcor!
         const Double_t  dt = kSmpRate * (             // from ns to samples
            ( (disCh-disXc) * kNgTopFern / kC_m_ns )   // from m to ns
            + dtcor );                                 // correct dt offset (ns)
         // FIXME: dt in samples, maxdt in ns
         if (TMath::Abs(dt)>maxdt) {
            // really don't like being out of bounds
            c = TMath::Exp(dt*dt);
            if (c>kReallyBig) {
               c = kReallyBig;
            }
         } else {
            // get the correlation coefficient for this dt (in num samples)
            c = gspc[ch][xc]->Eval(dt) - 1.0;
         }
         chi2 += c*c;
      }
   }
   return chi2;
}

Double_t getShiftLL(const Double_t* thetaPhi) {
   
   Double_t chi2(1), c(0), oo(0);

   TVector3 norm; norm.SetMagThetaPhi(1.0, thetaPhi[0], thetaPhi[1]);
   for (Int_t ch=0; ch<NSnConstants::kNchans; ++ch) {
      for (Int_t xc=0; xc<ch; ++xc) {
         Double_t dtcor=0;
         for (Int_t i=(ch-xc); i>0; --i) {
            dtcor += dtCorrs[ch-i];
         }
         const TVector3& posCh = getStnPos(ch);
         const TVector3& posXc = getStnPos(xc);
         const Double_t  disCh = -(posCh.Dot(norm));
         const Double_t  disXc = -(posXc.Dot(norm));
         // !!! check sign of delta(distance) and dtcor!
         Double_t  dt = 
            ( (disCh-disXc) * kNgTopFern / kC_m_ns ) // from m to ns
            + dtcor;                                 // correct dt offset (ns)
         const Double_t odt = dt;
         Bool_t oob=kFALSE;
         if (dt<-maxdt) {
            dt  = -maxdt;
            oob = kTRUE;
         } else if (dt>maxdt) {
            dt  = maxdt;
            oob = kTRUE;
         }
         dt *= kSmpRate;
         c = getProbFromCorrCoef(gspl[ch][xc]->Eval(dt));
         if (oob) {
            const Double_t wa = TMath::Abs(odt) - maxdt;
            oo += wa*wa;
         }
         /*
         if (TMath::Abs(dt)>maxdt) {
            // really don't like being out of bounds
            c = TMath::Exp(-dt*dt);
         } else {
            // get the correlation coefficient for this dt (in num samples)
            const Double_t  corco = gspl[ch][xc]->Eval(dt);
            c = getProbFromCorrCoef(corco);
         }
         if (c>1.0) {
            Fatal("getShiftLL","Got ll term > 1 (%g)",c);
         }
         */
         chi2 *= c;
      }
   }
   chi2 = -TMath::Log(chi2);
   chi2 += oo;
   return chi2;
}

TGraph GetSecondaryDecayCurve(int Points,const double & hi_ex_set,const double& decay_ex_set)
{
  TGraph curve;
  if(!gPrimaryReaction.IsSet()){
    std::cout<<"Reaction Masses have not been set"<<std::endl;
    exit(EXIT_FAILURE);
  }
  if(!gPrimaryReaction.BeamEnergy()){
    std::cout<<"Beam Energy has not been set"<<std::endl;
    exit(EXIT_FAILURE);
  }
 
  sim::RN_SimEvent evt1(gPrimaryReaction.BeamEnergy(),gPrimaryReaction.M_Beam(),gPrimaryReaction.M_Target(),gPrimaryReaction.M_Recoil(),gPrimaryReaction.M_Fragment());
  sim::RN_SimEvent evt2(gPrimaryReaction.M_Fragment(),gPrimaryReaction.M_Decay_Product(),gPrimaryReaction.M_Heavy_Decay());
  // Fill the points of the kinematic curve
  int p=0;
  while(p<Points){
    double theta_deg = 175; // assume backward angle from inverse kinematics
    double phi=2.*M_PI*global::myRnd.Rndm();
    
    TVector3 nv; nv.SetMagThetaPhi(1.,theta_deg*M_PI/180.0,phi);
    if(!evt1.radiate_in_CM(nv,hi_ex_set))
      continue;

    theta_deg = 180* p / Points;
    phi = 2.*M_PI*global::myRnd.Rndm();
    TVector3 pv; pv.SetMagThetaPhi(1.,theta_deg*M_PI/180.0,phi);
    if(!evt2.radiate_in_CM(evt1.getLVhi(),pv,decay_ex_set))
      continue;

    curve.SetPoint(p, evt2.getLVrad().Theta()*180/3.14,(double)(evt2.getLVrad().E() - evt2.getLVrad().M()));
    p++;
  }
  // end for(p)    

  return curve;
}

Double_t getAngleChi2(const Double_t* thetaPhi) {
   Double_t chi2(0), c(0);
   
   TVector3 norm; norm.SetMagThetaPhi(1.0, thetaPhi[0], thetaPhi[1]);
   for (Int_t ch=0; ch<NSnConstants::kNchans; ++ch) {
      for (Int_t xc=0; xc<ch; ++xc) {
         const TVector3& posCh = getStnPos(ch);
         const TVector3& posXc = getStnPos(xc);
         const Double_t  disCh = -(posCh.Dot(norm));
         const Double_t  disXc = -(posXc.Dot(norm));
         // !!! check sign of delta(distance) and dtcor!
         const Double_t  dt = (disCh-disXc) * kNgTopFern / kC_m_ns;  // from m to ns
         const Double_t bdt = thetaPhi[2+ch] - thetaPhi[2+xc];

         c = (dt-bdt);
         chi2 += c*c;
         /*
         Printf("(%g,%g) d[%d]=%g, d[%d]=%g, bdt=%g, sch=%g, sxc=%g, dt=%g, c=%g",
                thetaPhi[0]*TMath::RadToDeg(), thetaPhi[1]*TMath::RadToDeg(),
                ch, thetaPhi[2+ch], xc, thetaPhi[2+xc], bdt, disCh, disXc, dt, c);
         */
      }
   }
   
#ifdef DEBUG_ANGLE
   if (dbgth!=0) {
      dbgth->SetPoint(dbgth->GetN(), thetaPhi[0]*TMath::RadToDeg(),
                      chi2);
      dbgphi->SetPoint(dbgphi->GetN(), thetaPhi[1]*TMath::RadToDeg(),
                      chi2);
   }
#endif
   
   return chi2;
   
}

void fill(int const kf, TLorentzVector* b, double weight, TLorentzVector const& p1Mom, TLorentzVector const& p2Mom, TVector3 v00)
{
   int const centrality = floor(nCent * gRandom->Rndm());

   TVector3 const vertex = getVertex(centrality);
   // smear primary vertex
   // float const sigmaVertex = sigmaVertexCent[cent];
   // TVector3 const vertex(gRandom->Gaus(0, sigmaVertex), gRandom->Gaus(0, sigmaVertex), gRandom->Gaus(0, sigmaVertex));

   v00 += vertex;

   // smear momentum
   TLorentzVector const p1RMom = smearMom(0, p1Mom);
   TLorentzVector const p2RMom = smearMom(0, p2Mom);

   // smear position
   TVector3 const p1RPos = smearPosData(0, vertex.z(), centrality, p1RMom, v00);
   TVector3 const p2RPos = smearPosData(0, vertex.z(), centrality, p2RMom, v00);
   // TVector3 const kRPos = smearPos(kMom, kRMom, v00);
   // TVector3 const pRPos = smearPos(pMom, pRMom, v00);

   // reconstruct
   TLorentzVector const rMom = p1RMom + p2RMom;
   float const p1Dca = dca(p1Mom.Vect(), v00, vertex);
   float const p2Dca = dca(p2Mom.Vect(), v00, vertex);
   float const p1RDca = dca(p1RMom.Vect(), p1RPos, vertex);
   float const p2RDca = dca(p2RMom.Vect(), p2RPos, vertex);

   TVector3 v0;
   float const dca12 = dca1To2(p1RMom.Vect(), p1RPos, p2RMom.Vect(), p2RPos, v0);
   float const decayLength = (v0 - vertex).Mag();
   float const dcaD0ToPv = dca(rMom.Vect(), v0, vertex);
   float const cosTheta = (v0 - vertex).Unit().Dot(rMom.Vect().Unit());
   float const angle12 = p1RMom.Vect().Angle(p2RMom.Vect());

   TLorentzVector p1RMomRest = p1RMom;
   TVector3 beta;
   beta.SetMagThetaPhi(rMom.Beta(), rMom.Theta(), rMom.Phi());
   p1RMomRest.Boost(-beta);
   float const cosThetaStar = rMom.Vect().Unit().Dot(p1RMomRest.Vect().Unit());

   // save
   float arr[100];
   int iArr = 0;
   arr[iArr++] = centrality;
   arr[iArr++] = vertex.X();
   arr[iArr++] = vertex.Y();
   arr[iArr++] = vertex.Z();

   arr[iArr++] = kf;
   arr[iArr++] = b->M();
   arr[iArr++] = b->Perp();
   arr[iArr++] = b->PseudoRapidity();
   arr[iArr++] = b->Rapidity();
   arr[iArr++] = b->Phi();
   arr[iArr++] = v00.X();
   arr[iArr++] = v00.Y();
   arr[iArr++] = v00.Z();

   arr[iArr++] = rMom.M();
   arr[iArr++] = rMom.Perp();
   arr[iArr++] = rMom.PseudoRapidity();
   arr[iArr++] = rMom.Rapidity();
   arr[iArr++] = rMom.Phi();
   arr[iArr++] = v0.X();
   arr[iArr++] = v0.Y();
   arr[iArr++] = v0.Z();

   arr[iArr++] = dca12;
   arr[iArr++] = decayLength;
   arr[iArr++] = dcaD0ToPv;
   arr[iArr++] = cosTheta;
   arr[iArr++] = angle12;
   arr[iArr++] = cosThetaStar;

   arr[iArr++] = p1Mom.M();
   arr[iArr++] = p1Mom.Perp();
   arr[iArr++] = p1Mom.PseudoRapidity();
   arr[iArr++] = p1Mom.Rapidity();
   arr[iArr++] = p1Mom.Phi();
   arr[iArr++] = p1Dca;

   arr[iArr++] = p1RMom.M();
   arr[iArr++] = p1RMom.Perp();
   arr[iArr++] = p1RMom.PseudoRapidity();
   arr[iArr++] = p1RMom.Rapidity();
   arr[iArr++] = p1RMom.Phi();
   arr[iArr++] = p1RPos.X();
   arr[iArr++] = p1RPos.Y();
   arr[iArr++] = p1RPos.Z();
   arr[iArr++] = p1RDca;
   arr[iArr++] = tpcReconstructed(0,1,centrality,p1RMom);

   arr[iArr++] = p2Mom.M();
   arr[iArr++] = p2Mom.Perp();
   arr[iArr++] = p2Mom.PseudoRapidity();
   arr[iArr++] = p2Mom.Rapidity();
   arr[iArr++] = p2Mom.Phi();
   arr[iArr++] = p2Dca;

//.........这里部分代码省略.........

//.........这里部分代码省略.........

   Printf("generating events...");
   TStopwatch timer;
   timer.Start();
   
   for (UInt_t i=0; i<simevts; ++i) {
      
      if ( (i%1000)==0 ) {
         fprintf(stderr,"Processing %u/%u ...            \r",i,simevts);
      }
      
      // choose angles
      theta = (thetaOpt>360.) ? TMath::ACos( rnd.Uniform(-1.0, 0.0) ) 
                              : thetaOpt * TMath::DegToRad();
      phi   = (phiOpt>360.) ? rnd.Uniform(0.0, TMath::TwoPi())
                            : phiOpt * TMath::DegToRad();
      cone  = (coneOpt>360.) 
         ? rnd.Uniform(*(Cangs.begin()), *(Cangs.rbegin()))
         : coneOpt; // leave this one in degrees (as in the tree)
      CAng = findNearestAllowedAngle(Cangs, cone);
      
#ifdef DEBUG
      Printf("--- theta=%g, phi=%g, cone=%g",
             theta*TMath::RadToDeg(), phi*TMath::RadToDeg(), cone);
#endif
      
      // calculate channel shifts
      TArrayD pwdt = NSnChanCorl::GetPlaneWaveOffsets(theta,
                                                      phi,
                                                      zeros,
                                                      pos,
                                                      kNgTopFirn);
      TVector3 dir;
      dir.SetMagThetaPhi(1.0, theta, phi);
      
#ifdef DEBUG
      TObjArray graphs;
      graphs.SetOwner(kTRUE);
      TCanvas* c1 = new TCanvas("c1","c1",800,700);
      c1->Divide(2,2);
#endif
      
      for (UChar_t ch=0; ch<NSnConstants::kNchans; ++ch) {
         
         // look up the EAng, fhang for this antenna
         Float_t feang(0), fhang(0);
         findEangHang(nvec[ch], dir, feang, fhang);
         feang  = TMath::Abs(TVector2::Phi_mpi_pi(feang));
         fhang  = TMath::Abs(TVector2::Phi_mpi_pi(fhang));
         feang *= TMath::RadToDeg();
         fhang *= TMath::RadToDeg();
         // find closest allowed angle
         EAng[ch] = findNearestAllowedAngle(Eangs, feang);
         HAng[ch] = findNearestAllowedAngle(Hangs, fhang);
         const Long64_t ni = 
            nnt->GetEntryNumberWithIndex(EAng[ch] + (1000*HAng[ch]), CAng);
#ifdef DEBUG
         Printf("EAng=%g (%g), HAng=%g (%g), CAng=%g, ni=%lld",
                EAng[ch],feang,HAng[ch],fhang,CAng,ni);
#endif
         if (ni>-1) {
            nnt->GetEntry(ni);
#ifdef DEBUG
            c1->cd(ch+1);
            TGraph* och = wave->NewGraphForChan(0, kTRUE);
            const Int_t ochnp = och->GetN();

SimulationOutput* LightSimulator::Run(){
    
  TRandom3 randgen = TRandom3(0);

  if (debug) cout << ntoys*nrays << " light propagations to simulate" << endl;

  long counter_ray=0;

  if (output) delete output;
  output = new SimulationOutput(deposit);
    
  if (!isinsideboundaries()) return output;

  for (int nrun = 0; nrun<ntoys; nrun++){
            
    vector<Int_t> myphotons(4,0);

    for (int nray = 0; nray<nrays; nray++){

      counter_ray++;
      //	if (counter_ray%(ntoys*nrays/10)==0) cout << "Done " << counter_ray << " rays" << endl;

      Double_t phi = randgen.Uniform(-Pi(),Pi());
      Double_t costheta = randgen.Uniform(-1,1);
      TVector3 dir; dir.SetMagThetaPhi(1,ACos(costheta),phi);
      LightRay lr(deposit.position,dir);

      bool matched = false;
      Int_t matchx = 999;
      Int_t matchy = 999;
      Int_t firstmatchx = 999;
      Int_t firstmatchy = 999;

      Int_t index=0;
      MatchObject m;
      vector<MatchObject> matches;
      while (index>=0 && index<pars.max_distance_unit_propagation && index<Max(firstmatchx,firstmatchy)){
	if (propray(lr,kNS,index,matchx,m)){
	  if (!matched) firstmatchx=matchx;
	  matched=true;
	  matches.push_back(m);
	}
	if (propray(lr,kEW,index,matchy,m)){
	  if (!matched) firstmatchy=matchy;
	  matched=true;
	  matches.push_back(m);
	}
	index++;
      }

      vector<GoodMatchObject> goodmatches=convert_matches(matches);
	
      Double_t pathxy = pars.max_distance_unit_propagation*10*pars.xtalsize;
      GoodMatchObject finalmatch;
      for (size_t i=0; i<goodmatches.size(); i++){
	Double_t path = sqrt(pow(goodmatches[i].positionx-lr.origin.x(),2)+pow(goodmatches[i].positiony-lr.origin.y(),2));
	if (path<pathxy) {pathxy=path; finalmatch=goodmatches[i];}
      }
      if (pathxy>1e4){
	if (debug) cout << "LOOPER" << endl;
	continue;
      }

      Int_t channel = findchannel(finalmatch.x,finalmatch.y)-1;
     
      if (debug) cout << finalmatch.x << " " << finalmatch.y << " " << pathxy << " " << channel+1 << endl;

      Double_t path3d = pathxy/Sin(lr.dirvect.Theta());
     
      TVector3 rotated_origin_xy = lr.origin; rotated_origin_xy.SetZ(0);
      TVector3 rotated_impact = TVector3(finalmatch.positionx,finalmatch.positiony,0);
      TRotation r;
      r.RotateZ(-Pi()/4);
      rotated_origin_xy = r*rotated_origin_xy;
      rotated_impact = r*rotated_impact;

      Int_t nx = 0;
      Int_t ny = 0;
      Int_t nz = 0;
      {
	Int_t sx = finalmatch.x+finalmatch.y;
	if (sx==0) nx=0;
	else if (sx>0) nx = sx/2-1;
	else nx = TMath::Abs(sx)/2;
	Int_t dy = finalmatch.y-finalmatch.x;
	if (dy==0) ny=0;
	else if (dy>0) ny = dy/2-1;
	else ny = TMath::Abs(dy)/2;
	Float_t finalz = path3d*Cos(lr.dirvect.Theta())+deposit.position.z()-pars.xtalheight/2;
	nz = Int_t((fabs(finalz)+pars.xtalheight/2)/pars.xtalheight);
      }
      Int_t all_crossings = nx+ny+nz;

      Int_t my_paper_refl = 0;

      TVector2 difference = TVector2(fabs(rotated_origin_xy.x()-rotated_impact.x()),fabs(rotated_origin_xy.y()-rotated_impact.y()));
      if (difference.Phi()<limit_angle) my_paper_refl+=nx;
      if (Pi()/2-difference.Phi()<limit_angle) my_paper_refl+=ny;
      if (lr.dirvect.Theta()<limit_angle) my_paper_refl+=nz;

//.........这里部分代码省略.........

void RayTrace(AOpticsManager* manager, TCanvas* can3D)
{
  const int kNdeg = 8;
  TH2D* h2[kNdeg];
  TGraph* graph = new TGraph();
  TCanvas* can = new TCanvas("can", "can", 900, 900);
  TCanvas* can2= new TCanvas("can2", "can2", 900, 900);
  can->Divide(3, 3, 1e-10, 1e-10);

  TH2D* hMirror = new TH2D("hMirror", ";X (mm);Y (mm)", 1000, -7, 7, 1000, -7, 7);

  for(int i = 0; i < kNdeg; i++){

    double deg = i*0.5;
    TGeoTranslation raytr("raytr", -2*kF*TMath::Sin(deg*TMath::DegToRad()), 0, 2*kF*TMath::Cos(deg*TMath::DegToRad()));
    TVector3 dir;
    dir.SetMagThetaPhi(1, TMath::Pi() - deg*TMath::DegToRad(), 0);
    double lambda = 400*nm; // dummy
    ARayArray* array = ARayShooter::Square(lambda, 14*m, 401, 0, &raytr, &dir);

    manager->TraceNonSequential(*array);

    h2[i] = new TH2D("", Form("#it{#theta} = %3.1f#circ;x (mm); y (mm)", deg), 200, -40, 100, 200, -70, 70);
    TH2D tmp("", "", 100, -1e5, 1e5, 100, -1e5, 1e5);

    TObjArray* focused = array->GetFocused();

    for(Int_t j = 0; j <= focused->GetLast(); j++){
      ARay* ray = (ARay*)(*focused)[j];
      Double_t p[4];
      ray->GetLastPoint(p);
      tmp.Fill(p[0], p[1]);

      if (i == 0) {
        int n = ray->FindNodeNumberStartWith("mirror");
        const double* pn = ray->GetPoint(n);
        hMirror->Fill(pn[0]/m, pn[1]/m);
      } // if

      if (i == kNdeg - 1 && gRandom->Uniform() < 0.001) {
        TPolyLine3D* pol = ray->MakePolyLine3D();
        pol->SetLineColor(2);
        can3D->cd();
        pol->Draw();
      } // if
    } // j

    double meanx = tmp.GetMean();

    for(Int_t j = 0; j <= focused->GetLast(); j++){
      ARay* ray = (ARay*)(*focused)[j];
      Double_t p[4];
      ray->GetLastPoint(p);
      h2[i]->Fill((p[0] - meanx)/mm, p[1]/mm);
    } // j

    can->cd(i + 1);
    h2[i]->Draw("colz");

    if(i == 0){
      can2->cd();
      hMirror->Draw("colz");
    } // i

    delete array;
  } // i
}

//.........这里部分代码省略.........
            Double_t phi   = result[1];
            if (theta<0) {
               theta *= -1.0;
               phi   += TMath::Pi();
            }
            phi   = TVector2::Phi_0_2pi(phi);
            Printf("theta=%g, phi=%g (%g, %g)",
                   theta*TMath::RadToDeg(),
                   phi*TMath::RadToDeg(),
               result[0], result[1]);
            TCanvas* cdbgtp = new TCanvas("cdbgtp","cdbgtp",1000,1000);
            cdbgtp->Divide(2,3);
            Float_t xxsm[NSnConstants::kNsamps];
            for (Int_t i=0; i<NSnConstants::kNsamps; ++i) {
               xxsm[i] = i;
            }
            
            TGraph* gorg[NSnConstants::kNchans];
            TGraph* gocp[NSnConstants::kNchans];
            TGraph* gshf[NSnConstants::kNchans][NSnConstants::kNchans];
            for (Int_t ch=0; ch<NSnConstants::kNchans; ++ch) {
               gorg[ch] = new TGraph(NSnConstants::kNsamps, xxsm,
                                     &(filtered[ch][0]));
               gocp[ch] = new TGraph(NSnConstants::kNsamps, xxsm,
                                     &(filtered[ch][0]));
               for (Int_t xc=0; xc<ch; ++xc) {
                  Double_t dtcor=0;
                  for (Int_t i=(ch-xc); i>0; --i) {
                     dtcor += dtCorrs[ch-i];
                  }
                  const TVector3& posCh = getStnPos(ch);
                  const TVector3& posXc = getStnPos(xc);
                  TVector3 norm;
                  norm.SetMagThetaPhi(1.0, result[0], result[1]);
                  //norm.SetMagThetaPhi(1.0, 2.95833, 5.497787);
                  const Double_t  disCh = -(posCh.Dot(norm));
                  const Double_t  disXc = -(posXc.Dot(norm));
                  // !!! check sign of delta(distance) and dtcor!
                  const Double_t  dt = kSmpRate * (
                     ( (disCh-disXc) * kNgTopFern / kC_m_ns )   // from m to ns
                     + dtcor );
                  Printf("dt[%d,%d]=%g (dis[%d]=%g, dis[%d]=%g)",
                         ch,xc,dt, ch,disCh,xc,disXc);
                  gshf[ch][xc] = new TGraph(NSnConstants::kNsamps);
                  for (Int_t i=0; i<NSnConstants::kNsamps; ++i) {
                     gshf[ch][xc]->SetPoint(i, xxsm[i]+dt, filtered[xc][i]);
                  }
                  gshf[ch][xc]->SetMarkerStyle(7);
                  gshf[ch][xc]->SetLineColor(kRed);
                  gshf[ch][xc]->SetMarkerColor(kRed);
               }
               gorg[ch]->SetMarkerStyle(7);
               gorg[ch]->SetLineColor(kBlack);
               gorg[ch]->SetMarkerColor(kBlack);
               gocp[ch]->SetMarkerStyle(7);
               gocp[ch]->SetLineStyle(7);
               gocp[ch]->SetLineColor(kGreen+2);
               gocp[ch]->SetMarkerColor(kGreen+2);
            }
            gStyle->SetOptStat(0);
            Int_t pp=0;
            for (Int_t ch=0; ch<NSnConstants::kNchans; ++ch) {
               for (Int_t xc=0; xc<ch; ++xc) {
                  cdbgtp->cd(++pp);
                  hn = Form("hd_%d_%d",ch,xc);
                  TH2F* hd = new TH2F(hn.Data(),