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