Hi Eugene,
it is very hard to debug codes just from a plot, it is even harder to do
so on the new year's day. :)
all I can suggest is to have a stripped version of the code with only the
statements pertain to the generation of this plot.
it is bothersome that entry.twr is changes by codes other than the
GetEntry() line, if I read your message correctly, try protect the values
in entry that were read from the ntuple as strongly as possible and only
change the "derived" values such as entry.pe. for instance, avoid
SetAddress of any branches to &entry.pe.
if the stripped macro/code is less than 100 lines, then you can posted and
people will be more likely tolook at it after coming back from the
holidays. :)
I never saw any pattern like this.
Chi
On Sat, 31 Dec 2005, Eugene J. Geis wrote:
> Hi,
>
> I'm pulling out my hair. I'm using the flr files and putting together
> a truncated dataset using simple methods such as new TTree() and
> SetAddress()... blah blah. There is one point in my macro where I
> simply make a new variable for electron momentum and, depending
> on cuts, read in the value, recon.pe = entry.pwr, where entry.pwr
> is linked to the flr entry pwr, etc.
> I'm then doing kinematic calculations and writing in OTHER entries
> in my branches of recon and entry. Somehow, I'm getting a split in
> the data down the middle of my momentum distribution that
> looks like the included graph c1.ps where I've plotted pe:pp (p_elec vs.
> p_prot). This "software mistake" is actually rewritten all the way down
> into the pwr entry that is linked directly to the address of the flr pwr
> entry. When I look at the actual data in the flr files, they have no split like
> this. I'm completely lost... I've rewritten the macro three times and
> cleared out a load of code and I have no freaking clue. If anyone has a
> few minutes to see what the F&%$ I'm doing wrong, or what weird C++/
> ROOT error I might have haphazardly invoked, PLEASE let me know...
>
> eugene
>
> --------------------------------------------------------------------------
> Eugene Geis
> PhD Student, Physics Department, ASU
> Research Affiliate, MIT-Bates Laboratory of Nuclear Science
> eugene.geis@asu.edu
>
>
> const Float_t eBeam = 0.8500; //electron beam energy
>
> const Int_t nSects = 2 ; //number of sectors
> const Int_t nSpins = 6 ; //number of beam-target spin combinations
> const Int_t nX2Bins = 20 ;
>
> const Float_t mE = 0.000511; // mass of electron
> const Float_t mD = 1.875613; //mass of deuteron
> const Float_t mP = 0.938272; //mass of proton
> const Float_t mN = 0.939565; //mass of neutron
> const Float_t Vsc = 30.00000; // speed of light in vacuum
>
> const Float_t deg2rad = (4.0*atan(1.0))/180.0; //[rad] = deg2rad * [deg]
> const Float_t rad2deg = 180.0/(4.0*atan(1.0)); //[deg] = rad2deg * [rad]
>
> //*****************************************************************************
>
> struct Entry { //linked to flr entries
> Float_t ttl, ttr, atl, atr, ntl, ntr, ptl, ptr;
> Float_t pwl, twl, fwl, zwl, qwl, xwl;
> Float_t pwr, twr, fwr, zwr, qwr, xwr;
> Float_t ncl, ncr, acl, acr, tcl, tcr;
> Float_t Bwl, Bwr, trig;
> Float_t nrun;
> } entry; // variables inhereted from reconstruction
>
> struct Recon{
> Float_t Ee,pe,thetae,phie,betae;
> Float_t Ep,pp,thetap,phip,betap;
> Float_t pm,Em,Mm,Q2;
> Float_t pe_pp,pe_te,pe_tp;
> Float_t pp_pe,pp_tp,pp_te;
> Float_t te_pe,te_pp,te_tp;
> Float_t tp_pp,tp_pe,tp_te;
> Float_t Q2_pe,Q2_pp,Q2_te,Q2_tp,X2Q2,dQ2;
> Float_t pe_Q2,pp_Q2,te_Q2,tp_Q2;
> } recon; // calculated variables
>
> TTree *kc = new TTree("kc", "KineCorr NTuple"); // usefull for debuging
>
> Bool_t mc = 0; //true if the data is Monte Carlo
> Bool_t whgen = 0; //true if the data is white-generated
> Bool_t ECUT = 0; //true if I do nice momentum vs. theta cut on electron
> Int_t nEvents=999999999;
>
> Float_t omega,q,phiq,thetaq,Q2_center,tester,tester2;
> Int_t sect,npass;
> Float_t pnx,pny,pnz;
> Float_t qx,qy,qz;
> Float_t pmx,pmy,pmz;
> Float_t theta,phi;
> Float_t pB[4] = {0};
> Float_t pQ[4] = {0};
> Float_t test[6];
> Float_t Q2pe, Q2pp, Q2te, Q2tp, Q2center, d_Q2; //deleted function calling these
> Float_t emom,pmom,eth,pth; // one attempt at debugging.
> //*****************************************************************************
>
>
> TFile *fmc;
>
>
> //*****************************************************************************
>
> void InitAll(){
>
> gStyle->SetOptStat(0);
> gStyle->SetOptTitle(0);
> gStyle->SetOptDate(1);
> gStyle->SetOptFit(1);
>
> int arc=gApplication->Argc();
> char** arv=gApplication->Argv();
>
> // for(Int_t i=0;i<arc;i++){
> // TString arg(arv[i]);
> // if(arg=="-m"){mc=true;cout<<"Monte Carlo File"<<endl;}
> // else if((arg=="-n")&&(i+1<arc)){sscanf(arv[++i],"%i",&nEvents);}
> // else if (arg=="-white"){whgen=true;cout<<"White Generator: "<<whgen<<endl;}
> // else if (arg=="-cut"){ECUT=true;}
> // }
>
> gROOT->ProcessLine(".x flr.C");
>
> TMatrix ch = charge[0]; //array holding various charges
>
>
> if(ch(0,0)> 0){ // beam target
> printf("Charge in (beam-target) polarization states:\n");
> printf(" ( h, Pz, Pzz) charge\n");
> printf(" -------- -----------\n");
> printf(" ( +1, +1, +1) %10.4f\n",ch(0,0));
> printf(" ( -1, +1, +1) %10.4f\n",ch(1,0));
> printf(" ( +1, -1, +1) %10.4f\n",ch(0,1));
> printf(" ( -1, -1, +1) %10.4f\n",ch(1,1));
> printf(" ( +1, 0, -2) %10.4f\n",ch(0,2));
> printf(" ( -1, 0, -2) %10.4f\n",ch(1,2));
> //printf(" ( -0, +0) %10.4f %10.4f\n",ch(0,2), ch(1,2));
> printf("\n");
> }
> else{
> printf("\nERROR: NO CHARGE\n\n");
> exit(1);
> }
>
>
> SetAddress(chain,&entry);
> kc->Branch("entry",&entry,"ttl:ttr:atl:atr:ntl:ntr:ptl:ptr:" // set TTree entries
> "pwl:twl:fwl:zwl:qwl:xwl:"
> "pwr:twr:fwr:zwr:qwr:xwr:"
> "ncl:ncr:acl:acr:tcl:tcr:"
> "Bwl:Bwr:trig:nrun");
> kc->Branch("recon",&recon,"Ee:pe:thetae:phie:betae:"
> "Ep:pp:thetap:phip:betap:pm:Em:Mm:Q2:"
> "pe_pp:pe_te:pe_tp:"
> "pp_pe:pp_tp:pp_te:"
> "te_pe:te_pp:te_tp:"
> "tp_pp:tp_pe:tp_te:"
> "Q2_pe:Q2_pp:Q2_te:Q2_tp:X2Q2:dQ2:"
> "pe_Q2:pp_Q2:te_Q2:tp_Q2");
>
> }
>
> int SetAddress(TChain** chain,struct Entry* entry) { // link to flr
> chain[0]->SetBranchStatus("ttl", 1);
> chain[0]->SetBranchStatus("ttr", 1);
> chain[0]->SetBranchStatus("atl", 1);
> chain[0]->SetBranchStatus("atr", 1);
> chain[0]->SetBranchStatus("ntl", 1);
> chain[0]->SetBranchStatus("ntr", 1);
> chain[0]->SetBranchStatus("ptl", 1);
> chain[0]->SetBranchStatus("ptr", 1);
>
> chain[0]->SetBranchStatus("pwl", 1);
> chain[0]->SetBranchStatus("twl", 1);
> chain[0]->SetBranchStatus("fwl", 1);
> chain[0]->SetBranchStatus("zwl", 1);
> chain[0]->SetBranchStatus("qwl", 1);
> chain[0]->SetBranchStatus("xwl", 1);
>
> chain[0]->SetBranchStatus("pwr", 1);
> chain[0]->SetBranchStatus("twr", 1);
> chain[0]->SetBranchStatus("fwr", 1);
> chain[0]->SetBranchStatus("zwr", 1);
> chain[0]->SetBranchStatus("qwr", 1);
> chain[0]->SetBranchStatus("xwr", 1);
>
> chain[0]->SetBranchStatus("ncl", 1);
> chain[0]->SetBranchStatus("ncr", 1);
> chain[0]->SetBranchStatus("acl", 1);
> chain[0]->SetBranchStatus("acr", 1);
> chain[0]->SetBranchStatus("tcl", 1);
> chain[0]->SetBranchStatus("tcr", 1);
>
> chain[0]->SetBranchStatus("Bwl", 1);
> chain[0]->SetBranchStatus("Bwr", 1);
>
> chain[0]->SetBranchStatus("trig", 1);
> chain[0]->SetBranchStatus("nrun", 1);
>
> chain[0]->SetBranchAddress("ttl", &entry->ttl);
> chain[0]->SetBranchAddress("ttr", &entry->ttr);
> chain[0]->SetBranchAddress("atl", &entry->atl);
> chain[0]->SetBranchAddress("atr", &entry->atr);
> chain[0]->SetBranchAddress("ntl", &entry->ntl);
> chain[0]->SetBranchAddress("ntr", &entry->ntr);
> chain[0]->SetBranchAddress("ptl", &entry->ptl);
> chain[0]->SetBranchAddress("ptr", &entry->ptr);
>
> chain[0]->SetBranchAddress("pwl", &entry->pwl);
> chain[0]->SetBranchAddress("twl", &entry->twl);
> chain[0]->SetBranchAddress("fwl", &entry->fwl);
> chain[0]->SetBranchAddress("zwl", &entry->zwl);
> chain[0]->SetBranchAddress("qwl", &entry->qwl);
> chain[0]->SetBranchAddress("xwl", &entry->xwl);
>
> chain[0]->SetBranchAddress("pwr", &entry->pwr);
> chain[0]->SetBranchAddress("twr", &entry->twr);
> chain[0]->SetBranchAddress("fwr", &entry->fwr);
> chain[0]->SetBranchAddress("zwr", &entry->zwr);
> chain[0]->SetBranchAddress("qwr", &entry->qwr);
> chain[0]->SetBranchAddress("xwr", &entry->xwr);
>
> chain[0]->SetBranchAddress("ncl", &entry->ncl);
> chain[0]->SetBranchAddress("ncr", &entry->ncr);
> chain[0]->SetBranchAddress("acl", &entry->acl);
> chain[0]->SetBranchAddress("acr", &entry->acr);
> chain[0]->SetBranchAddress("tcl", &entry->tcl);
> chain[0]->SetBranchAddress("tcr", &entry->tcr);
>
> chain[0]->SetBranchAddress("Bwl", &entry->Bwl);
> chain[0]->SetBranchAddress("Bwr", &entry->Bwr);
>
> chain[0]->SetBranchAddress("trig", &entry->trig);
> chain[0]->SetBranchAddress("nrun", &entry->nrun);
>
> }
>
> //*****************************************************************************
>
> bool LooseCut() {
> return (fabs(entry.zwl) < 20 && fabs(entry.zwr) < 20 && entry.trig==1);
> }
>
> bool eLeftpRight() {
> return (entry.qwl == -1 && entry.ntl >-1 && entry.ntr > -1 && entry.qwr == 1);
> }
>
> bool pLefteRight() {
> return (entry.qwr == -1 && entry.ntr >-1 && entry.ntl > -1 && entry.qwl == 1);
> }
>
> bool MmCut(Float_t tol){
> return (fabs(recon.Mm) < tol);
> }
>
> bool CopCut(Float_t tol){
> return (fabs(entry.fwl-entry.fwr+180.0) < tol);
> }
>
> Bool_t BetaCut(Float_t beta,Float_t min,Float_t max){
> return((beta>min)&&(beta<max));
> }
>
> //*****************************************************************************
>
> void Blast2Q(Float_t *pB, Float_t *pQ, Float_t theta, Float_t phi) { // rotation
>
> pQ[0] = pB[0]*(-cos(theta*deg2rad)*cos(phi*deg2rad))
> + pB[1]*(-cos(theta*deg2rad)*sin(phi*deg2rad))
> + pB[2]*sin(theta*deg2rad);
> pQ[1] = pB[0]*sin(phi*deg2rad) + pB[1]*(-cos(phi*deg2rad));
> pQ[2] = pB[0]*(sin(theta*deg2rad)*cos(phi*deg2rad))
> + pB[1]*(sin(theta*deg2rad)*sin(phi*deg2rad))
> + pB[2]*cos(theta*deg2rad);
> pQ[3] = pB[3];
> }
>
> //**************************************************************************
> void KineCalc() { /// calculate kinematics
> recon.Ee = sqrt(emom*emom+mE*mE);
>
> recon.Ep = sqrt(pmom*pmom+mP*mP);
>
> pB[0] = pnx = pmom*sin(pth*deg2rad)*cos(recon.phip*deg2rad);
> pB[1] = pny = pmom*sin(pth*deg2rad)*sin(recon.phip*deg2rad);
> pB[2] = pnz = pmom*cos(pth*deg2rad);
>
> qx = -emom*sin(eth*deg2rad)*cos(recon.phie*deg2rad);
> qy = -emom*sin(eth*deg2rad)*sin(recon.phie*deg2rad);
> qz = eBeam - emom*cos(eth*deg2rad);
> omega = eBeam-recon.Ee;
>
> q = sqrt(qx*qx+qy*qy+qz*qz);
> phiq = (abs(recon.phie) < 90) ? 180 + recon.phie : recon.phie - 180;
> thetaq = acos(qz/q)*rad2deg;
>
> Blast2Q(pB,pQ,thetaq,phiq);
>
> pmx = -pQ[0];
> pmy = -pQ[1];
> pmz = q - pQ[2];
> recon.pm = sqrt(q*q - 2*q*pQ[2] + pmom*pmom);
> recon.Em = (mP+eBeam-recon.Ee-recon.Ep);
> recon.Mm = sqrt(recon.Em*recon.Em-recon.pm*recon.pm);
> recon.Q2 = q*q-omega*omega;
>
> if (fabs(recon.Mm)<0.1) { // calculate elastic relations
>
> recon.pe_pp = mP+eBeam-sqrt(pmom**2+mP**2);
>
> recon.pe_te = (eBeam*mP+mE**2)/(mP+eBeam-eBeam*cos(eth*deg2rad));
>
> recon.pp_pe = sqrt((mP+eBeam-emom)**2-mP**2);
>
> recon.pp_tp = 2*mP*(eBeam+mP)/((eBeam*cos(pth*deg2rad))*(((eBeam+mP)/
> (eBeam*cos(pth*deg2rad)))**2-1));
>
> recon.pe_tp = mP+eBeam-sqrt(recon.pp_tp**2+mP**2);
>
> recon.pp_te = sqrt((mP+eBeam-recon.pe_te)**2-mP**2);
>
> test[0] = (emom*(eBeam+mP)-eBeam*mP-mE**2)/(eBeam*emom);
> test[1] = ((mP+eBeam-sqrt(pmom**2+mP**2))*(eBeam+mP)-eBeam*mP-mE**2)/
> (eBeam*(mP+eBeam-sqrt(pmom**2+mP**2)));
> test[2] = (recon.pe_tp*(eBeam+mP)-eBeam*mP-mE**2)/(eBeam*recon.pe_tp);
> test[3] = (eBeam+mP)*(sqrt(pmom**2+mP**2)-mP)/(eBeam*pmom);
> test[4] = (eBeam+mP)*(sqrt((sqrt((mP+eBeam-emom)**2-mP**2))**2+mP**2)-mP)/
> (eBeam*(sqrt((mP+eBeam-emom)**2-mP**2)));
> test[5] = (eBeam+mP)*(sqrt(recon.pp_te**2+mP**2)-mP)/(eBeam*recon.pp_te);
>
>
>
> if (fabs(test[0])<=1.0) recon.te_pe = rad2deg*acos(test[0]); // if acos() can be invoked
> else recon.te_pe = -1;
>
> if (fabs(test[1])<=1.0) recon.te_pp = rad2deg*acos(test[1]);
> else recon.te_pp = -1;
>
> if (fabs(test[2])<=1.0) recon.te_tp = rad2deg*acos(test[2]);
> else recon.te_tp = -1;
>
> if (fabs(test[3])<=1.0) recon.tp_pp = rad2deg*acos(test[3]);
> else recon.tp_pp = -1;
>
> if (fabs(test[4])<=1.0) recon.tp_pe = rad2deg*acos(test[4]);
> else recon.tp_pe = -1;
>
> if (fabs(test[5])<=1.0) recon.tp_te = rad2deg*acos(test[5]);
> else recon.tp_te = -1;
>
> if (recon.te_pe>0&&recon.te_pp>-1&&recon.te_tp>-1) { // if reasonable values are found
> recon.Q2_pe = 4*0.85*emom*(sin(recon.te_pe*deg2rad/2.0))**2;
> recon.Q2_pp = 4*0.85*recon.pe_pp*(sin(recon.te_pp*deg2rad/2.0))**2;
> recon.Q2_te = 4*0.85*recon.pe_te*(sin(eth*deg2rad/2.0))**2;
> recon.Q2_tp = 4*0.85*recon.pe_tp*(sin(recon.te_tp*deg2rad/2.0))**2;
> }
> kc->Fill(); // fill ntuple
> npass++;
> }
> }
>
>
>
> void show_corr() {
>
> InitAll(); // Initialize ntuples and data
> npass = 0;
> for(Int_t i=0;((i<chain[0]->GetEntries())&&(i<nEvents));++i){
> if(i%50000==0){printf("\n%9d/%6d/ ",i,npass);}
> else if(i%10000==0) cout<<":"<<flush;
> else if(i%1000==0) cout<<"."<<flush;
>
> sect = -1;
> chain[0]->GetEntry(i);
>
> if ( LooseCut() ) { // vertex
> if ( eLeftpRight() ){ // sector dependent ep-elastic
>
> emom = recon.pe = entry.pwl; // MOMENTUM with global variable attempt at debug
> eth = recon.thetae = entry.twl;
> recon.phie = entry.fwl;
> recon.betae = entry.Bwl;
> pmom = recon.pp = entry.pwr; // MOMENTUM
> pth = recon.thetap = entry.twr;
> recon.phip = entry.fwr;
> recon.betap = entry.Bwr;
>
> sect = 0;
> }
> if ( pLefteRight() ){
> emom = recon.pe = entry.pwr; //MOMENTUM
> eth = recon.thetae = entry.twr;
> recon.phie = entry.fwr;
> recon.betae = entry.Bwr;
> pmom = recon.pp = entry.pwl; //MOMENTUM
> pth = recon.thetap = entry.twl;
> recon.phip = entry.fwl;
> recon.betap = entry.Bwl;
>
> sect = 1;
>
> }
>
> if(sect>-1&&CopCut(5.0)&&recon.betae>0.99){ // IF passed cuts, calculate
> KineCalc(); // calculate kinematics
> }
>
> }
> }
> cout << "I'm calling output file."<<endl;
>
> fout = new TFile("/scratch/bud08/egeis/data/kc-1.root","recreate"); // WRITE FILE.
> cout << "I'm writing kc."<<endl;
> kc->Write();
> cout << "I'm closing the output file."<<endl;
> fout->Close();
>
> }
>
>
This archive was generated by hypermail 2.1.2 : Mon Feb 24 2014 - 14:07:32 EST