Thanks everybody for taking a look... I don't completely understand it but I think there's a problem in
the way I'm calculating missing mass and then taking the square root is creating some kind of slice in
my momenta, probably at the point where missing mass perfectly equals zero; Mm = 0. I've been
using the MakeClass() method for several things and I have written my macro both ways with equal
results. (just finished crunching another....) It's definitely the missing mass calculation. Everything is
cured!! woohoo. Dios Mio. I haven't had the problem with my neutron data, but maybe I should take a
look because I think I calculate missing mass the same way that created my error in this program.
Granted, the proton mass >> 0 probably saves any asymptotic error in the missing mass calculation in
the e'n kinematics, but who knows... better safe than sorry. Better, updated kinematic corrections
coming to an email near you...
thanks again.
-eugene
Quoting Christopher Crawford <chris2@lns.mit.edu>:
> Hi Eugene,
> On the order of branches, I wouldn't even attempt to create my own
>
> ntuple event structure; it is much easier and less error-prone to
> just load the ntuple 'flr' or whatever, and then issue the command
>
> flr->MakeClass("flr_t");
>
> Then you just load "flr_t.h", and use the methods of that class to
> loop through the events. That way if the class structure changes,
> you just reissue the command, and your code still works. There is an
>
> example class in "blast/exp/analysis/utils/lr_t.h".
> --Chris
>
>
> On Jan 1, 2006, at 13:58:09, Nikolas Meitanis wrote:
>
> > Hey Eugene,
> > I looked at the macro (didnt run it) and have to admit I dont see
> > anything wrong with it...
> > Have you changed anything recently in this , assuming you have run
> > it before?
> >
> > Have you tried piping in events from the flr to your ntuple without
>
> > cuts etc?
> > I suppose there is a problem either with the ntuple construction,
> > the Struct , SetBranch etc,
> > or with the cuts. I didnt see anything wrong with the former. I
> > have been using a similar macro
> > and the only problems I have ever had were related to messing up
> > the order of variables
> > (I think they have to be in the same order in the
> >
> > 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
> >
> >
> > AND in the Branch:
> >
> > 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");
> >
> > }
> >
> > , and in this case they seem ok).
> > Another problem was related to the length of the Branch (somehow it
>
> > wouldnt write more than
> > a certain number of variables, have you tried using fewer variables
>
> > in your recon? There may be
> > a limit...).
> >
> > Other than that, the problem could be in the cuts:
> >
> > if ( LooseCut() ) {
> > if ( eLeftpRight() ){ ...
> > }
> > if ( pLefteRight() ){
> > ...
> > }
> >
> > }
> > etc. I am a little suspicious of the double equality thing, which I
>
> > have used before too but seemed unreliable somehow...
> >
> > emom = recon.pe = entry.pwr;
> > I would try without that too and see what happens.
> >
> > I cant see anything obvious though.
> >
> > n
> >
> >
> >
> >
> >
> >
> >
> >
> >
> > 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");
> >>
> >
> >
> > 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();
> >>
> >> }
> >>
> >>
>
>
--------------------------------------------------------------------------
Eugene Geis
PhD Student, Physics Department, ASU
Research Affiliate, MIT-Bates Laboratory of Nuclear Science
eugene.geis@asu.edu
--------------------------------------------------------------------------
http://quickreaction.blogspot.com
This archive was generated by hypermail 2.1.2 : Mon Feb 24 2014 - 14:07:32 EST