Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-05 08:12:06

0001 #include <vector>
0002 int eicsmear_dvmp_tree(TString filename_mc = "/sphenix/user/gregtom3/data/Fall2018/JPsi_reco_studies/reconstructedQ2/18x275_DVMP_1M_ascii_converted.root",
0003                TString filename_mc_smeared = "/sphenix/user/gregtom3/data/Fall2018/JPsi_reco_studies/reconstructedQ2/18x275_DVMP_1M_ascii_converted.smear.root")
0004 {
0005 
0006   /* PRELIMINARY ROOT STUFF */
0007 
0008   /* Loading libraries and setting sphenix style */
0009   gSystem->Load("libeicsmear");
0010   gROOT->LoadMacro("./sPHENIXStyle/sPhenixStyle.C");
0011   SetsPhenixStyle();
0012 
0013 
0014 
0015   /* INPUT */
0016 
0017 
0018 
0019   /* Input files */
0020   TFile *file_mc = new TFile(filename_mc, "OPEN");
0021   TFile *file_mc_smeared = new TFile(filename_mc_smeared, "OPEN");
0022   
0023   /* Trees */
0024   TTree *tree = (TTree*)file_mc->Get("EICTree");
0025   Double32_t Q2fromEICTree;
0026   tree->SetBranchAddress("trueQ2",&Q2fromEICTree);
0027   TTree *tree_smeared= (TTree*)file_mc_smeared->Get("Smeared");
0028 
0029   /* Output File */
0030   TFile *myfile = new TFile("/sphenix/user/gregtom3/data/Fall2018/JPsi_reco_studies/reconstructedQ2/analysisTree.root","RECREATE");
0031   /* Add friend */
0032   tree->AddFriend(tree_smeared);
0033 
0034 
0035   /* Create data tree because this is getting hectic */
0036   TTree *theTree = new TTree("DVMP_Tree","");
0037   Float_t de_phi_truth, de_eta_truth, de_pt_truth, de_p_truth;
0038   Float_t dp_phi_truth, dp_eta_truth, dp_pt_truth, dp_p_truth;
0039   Float_t sp_phi_truth, sp_eta_truth, sp_pt_truth, sp_p_truth;
0040   Float_t se_phi_truth, se_eta_truth, se_pt_truth, se_p_truth;
0041   Float_t jpsi_phi_truth, jpsi_eta_truth, jpsi_pt_truth, jpsi_p_truth;
0042   Float_t de_phi_reco, de_eta_reco, de_pt_reco, de_p_reco;
0043   Float_t dp_phi_reco, dp_eta_reco, dp_pt_reco, dp_p_reco;
0044   Float_t sp_phi_reco, sp_eta_reco, sp_pt_reco, sp_p_reco;
0045   Float_t se_phi_reco, se_eta_reco, se_pt_reco, se_p_reco;
0046   Float_t jpsi_phi_reco, jpsi_eta_reco, jpsi_pt_reco, jpsi_p_reco;
0047   Float_t invariant_reco_decay;
0048   Float_t invariant_reco_scatter;
0049   Float_t jpsi_px_truth, jpsi_py_truth, jpsi_pz_truth;
0050   Float_t jpsi_px_reco, jpsi_py_reco, jpsi_pz_reco;
0051   Double_t Q2;
0052   
0053   Float_t t;
0054   theTree->Branch("Q2",&Q2,"Event QSquared");
0055   theTree->Branch("t",&t,"Event t");
0056   theTree->Branch("de_phi_truth",&de_phi_truth,"Decay Electron Truth Phi");
0057   theTree->Branch("de_eta_truth",&de_eta_truth,"Decay Electron Truth Eta");
0058   theTree->Branch("de_pt_truth",&de_pt_truth,"Decay Electron Truth Pt");
0059   theTree->Branch("de_p_truth",&de_p_truth,"Decay Electron Truth P");
0060   theTree->Branch("dp_phi_truth",&dp_phi_truth,"Decay Positron Truth Phi");
0061   theTree->Branch("dp_eta_truth",&dp_eta_truth,"Decay Positron Truth Eta");
0062   theTree->Branch("dp_pt_truth",&dp_pt_truth,"Decay Positron Truth Pt");
0063   theTree->Branch("dp_p_truth",&dp_p_truth,"Decay Positron Truth P");
0064   theTree->Branch("sp_phi_truth",&sp_phi_truth,"Scattered Proton Truth Phi");
0065   theTree->Branch("sp_eta_truth",&sp_eta_truth,"Scattered Proton Truth Eta");
0066   theTree->Branch("sp_pt_truth",&sp_pt_truth,"Scattered Proton Truth Pt");
0067   theTree->Branch("sp_p_truth",&sp_p_truth,"Scattered Proton Truth P");
0068   theTree->Branch("se_phi_truth",&se_phi_truth,"Scattered Electron Truth Phi");
0069   theTree->Branch("se_eta_truth",&se_eta_truth,"Scattered Electron Truth Eta");
0070   theTree->Branch("se_pt_truth",&se_pt_truth,"Scattered Electron Truth Pt");
0071   theTree->Branch("se_p_truth",&se_p_truth,"Scattered Electron Truth P");
0072   theTree->Branch("jpsi_phi_truth",&jpsi_phi_truth,"Jpsi Truth Phi");
0073   theTree->Branch("jpsi_eta_truth",&jpsi_eta_truth,"Jpsi Truth Eta");
0074   theTree->Branch("jpsi_pt_truth",&jpsi_pt_truth,"Jpsi Truth Pt");
0075   theTree->Branch("jpsi_p_truth",&jpsi_p_truth,"Jpsi Truth P");
0076   theTree->Branch("de_phi_reco",&de_phi_reco,"Decay Electron Reco Phi");
0077   theTree->Branch("de_eta_reco",&de_eta_reco,"Decay Electron Reco Eta");
0078   theTree->Branch("de_pt_reco",&de_pt_reco,"Decay Electron Reco Pt");
0079   theTree->Branch("de_p_reco",&de_p_reco,"Decay Electron Reco P");
0080   theTree->Branch("dp_phi_reco",&dp_phi_reco,"Decay Positron Reco Phi");
0081   theTree->Branch("dp_eta_reco",&dp_eta_reco,"Decay Positron Reco Eta");
0082   theTree->Branch("dp_pt_reco",&dp_pt_reco,"Decay Positron Reco Pt");
0083   theTree->Branch("dp_p_reco",&dp_p_reco,"Decay Positron Reco P");
0084   theTree->Branch("sp_phi_reco",&sp_phi_reco,"Scattered Proton Reco Phi");
0085   theTree->Branch("sp_eta_reco",&sp_eta_reco,"Scattered Proton Reco Eta");
0086   theTree->Branch("sp_pt_reco",&sp_pt_reco,"Scattered Proton Reco Pt");
0087   theTree->Branch("sp_p_reco",&sp_p_reco,"Scattered Proton Reco P");
0088   theTree->Branch("se_phi_reco",&se_phi_reco,"Scattered Electron Reco Phi");
0089   theTree->Branch("se_eta_reco",&se_eta_reco,"Scattered Electron Reco Eta");
0090   theTree->Branch("se_pt_reco",&se_pt_reco,"Scattered Electron Reco Pt");
0091   theTree->Branch("se_p_reco",&se_p_reco,"Scattered Electron Reco P");
0092   theTree->Branch("jpsi_phi_reco",&jpsi_phi_reco,"Jpsi Reco Phi");
0093   theTree->Branch("jpsi_eta_reco",&jpsi_eta_reco,"Jpsi Reco Eta");
0094   theTree->Branch("jpsi_pt_reco",&jpsi_pt_reco,"Jpsi Reco Pt");
0095   theTree->Branch("jpsi_p_reco",&jpsi_p_reco,"Jpsi Reco P");
0096   theTree->Branch("jpsi_px_truth",&jpsi_px_truth,"Jpsi Truth Px");
0097   theTree->Branch("jpsi_py_truth",&jpsi_py_truth,"Jpsi Truth Py");
0098   theTree->Branch("jpsi_pz_truth",&jpsi_pz_truth,"Jpsi Truth Pz");
0099   theTree->Branch("jpsi_px_reco",&jpsi_px_reco,"Jpsi Reco Px");
0100   theTree->Branch("jpsi_py_reco",&jpsi_py_reco,"Jpsi Reco Py");
0101   theTree->Branch("jpsi_pz_reco",&jpsi_pz_reco,"Jpsi Reco Pz");
0102   theTree->Branch("invariant_reco_decay",&invariant_reco_decay,"Invariant Mass Decay Particles");
0103   theTree->Branch("invariant_reco_scatter",&invariant_reco_scatter,"Invariant Mass Scattered Electron");
0104   
0105   /* Momentum cut for plotting */
0106   TCut pcut("Smeared.particles.p<=50&&Smeared.particles.p>1");
0107   
0108   /* Eta cut for plotting */
0109   TCut ecut("-log(tan(Smeared.particles.theta/2))<=4&&-log(tan(Smeared.particles.theta/2))>=-4");
0110 
0111   /* Scattered lepton cut */
0112   TCut cut_scattered_lepton("particles.id==11&&particles.KS==2&&Smeared.particles.p<=50&&Smeared.particles.p>1");
0113 
0114   /* Decay electron cut */
0115   TCut cut_decay_electron("particles.id==11&&particles.KS==1&&Smeared.particles.p<=50&&Smeared.particles.p>1");
0116 
0117   /* Decay positron cut */
0118   TCut cut_decay_positron("particles.id==-11&&particles.KS==1&&Smeared.particles.p<=50&&Smeared.particles.p>1");
0119 
0120 
0121 
0122   /* PLOTTING */
0123 
0124   erhic::EventMilou *event  = NULL;
0125   Smear::Event       *eventS = NULL;
0126 
0127   tree->SetBranchAddress("event",&event);
0128   tree_smeared->SetBranchAddress("eventS", &eventS);
0129   
0130   unsigned max_event = tree_smeared->GetEntries();
0131   std::vector<float> particle_eta;
0132       particle_eta.push_back(0);
0133       particle_eta.push_back(0);
0134       particle_eta.push_back(0);
0135       particle_eta.push_back(0);
0136       particle_eta.push_back(0);
0137       std::vector<float> particle_theta;
0138       particle_theta.push_back(0);
0139       particle_theta.push_back(0);
0140       particle_theta.push_back(0);
0141       particle_theta.push_back(0);
0142       particle_theta.push_back(0);
0143       std::vector<float> particle_phi;
0144       particle_phi.push_back(0);
0145       particle_phi.push_back(0);
0146       particle_phi.push_back(0);
0147       particle_phi.push_back(0);
0148       particle_phi.push_back(0);
0149       std::vector<float> particle_pt;
0150       particle_pt.push_back(0);
0151       particle_pt.push_back(0);
0152       particle_pt.push_back(0);
0153       particle_pt.push_back(0);
0154       particle_pt.push_back(0);
0155       std::vector<float> particle_p;
0156       particle_p.push_back(0);
0157       particle_p.push_back(0);
0158       particle_p.push_back(0);
0159       particle_p.push_back(0);
0160       particle_p.push_back(0);
0161       
0162       std::vector<float> true_particle_eta;
0163       true_particle_eta.push_back(0);
0164       true_particle_eta.push_back(0);
0165       true_particle_eta.push_back(0);
0166       true_particle_eta.push_back(0);
0167       true_particle_eta.push_back(0);
0168       std::vector<float> true_particle_theta;
0169       true_particle_theta.push_back(0);
0170       true_particle_theta.push_back(0);
0171       true_particle_theta.push_back(0);
0172       true_particle_theta.push_back(0);
0173       true_particle_theta.push_back(0);
0174       std::vector<float> true_particle_phi;
0175       true_particle_phi.push_back(0);
0176       true_particle_phi.push_back(0);
0177       true_particle_phi.push_back(0);
0178       true_particle_phi.push_back(0);
0179       true_particle_phi.push_back(0);
0180       std::vector<float> true_particle_p;
0181       true_particle_p.push_back(0);
0182       true_particle_p.push_back(0);
0183       true_particle_p.push_back(0);
0184       true_particle_p.push_back(0);
0185       true_particle_p.push_back(0);
0186       std::vector<float> true_particle_pt;
0187       true_particle_pt.push_back(0);
0188       true_particle_pt.push_back(0);
0189       true_particle_pt.push_back(0);
0190       true_particle_pt.push_back(0);
0191       true_particle_pt.push_back(0);
0192   for ( unsigned ievent = 0; ievent < max_event/1000; ievent++ )
0193     {
0194       if ( ievent%1000 == 0 )
0195         cout << "Processing event " << ievent << endl;
0196 
0197       /* load event */
0198       tree->GetEntry(ievent);
0199       tree_smeared->GetEntry(ievent);
0200       cout << Q2fromEICTree << endl;
0201      
0202       unsigned ntracks = eventS->GetNTracks();
0203       
0204       
0205 
0206       
0207       for ( unsigned itrack = 0; itrack < ntracks; itrack++ )
0208         {
0209       
0210       Smear::ParticleMCS * sparticle = eventS->GetTrack( itrack );
0211       erhic::ParticleMC * tparticle = event->GetTrack(itrack);
0212       int type_of_particle = -1;
0213 
0214       // Decay Electron
0215       if(tparticle->Id().Code()==11&&tparticle->GetStatus()==1)
0216         {
0217           type_of_particle = 0;
0218         }
0219 
0220       // Scattered Electron
0221       if(tparticle->Id().Code()==11&&tparticle->GetStatus()==2)
0222         {
0223           type_of_particle = 1;
0224         }
0225       
0226       // Decay Positron
0227       if(tparticle->Id().Code()==-11&&tparticle->GetStatus()==1)
0228         {
0229           type_of_particle = 2;
0230         }
0231 
0232       // Scattered Proton
0233       if(tparticle->Id().Code()==2212&&tparticle->GetStatus()==1)
0234         {
0235           type_of_particle = 3;
0236         }
0237       // JPsi
0238       if(tparticle->Id().Code()==443)
0239         {
0240           type_of_particle = 4;
0241         }
0242       
0243       // Test that the smear particle was found
0244       if(sparticle==NULL&&type_of_particle!=4)
0245         {
0246           particle_eta.at(type_of_particle)=NULL;
0247           particle_theta.at(type_of_particle)=NULL;
0248           particle_phi.at(type_of_particle)=NULL;
0249           particle_pt.at(type_of_particle)=NULL;
0250           particle_p.at(type_of_particle)=NULL;
0251           
0252           true_particle_eta.at(type_of_particle)=tparticle->GetEta();
0253           true_particle_theta.at(type_of_particle)=tparticle->GetTheta();
0254           true_particle_phi.at(type_of_particle)=tparticle->GetPhi();
0255           true_particle_p.at(type_of_particle)=tparticle->GetP();
0256           true_particle_pt.at(type_of_particle)=tparticle->GetPt();
0257 
0258           //true_particle_eta.at(type_of_particle)=NULL;
0259           //true_particle_phi.at(type_of_particle)=NULL;
0260           //true_particle_p.at(type_of_particle)=NULL;
0261           //true_particle_pt.at(type_of_particle)=NULL;
0262         }
0263       else if(type_of_particle!=4)
0264         {
0265           
0266           particle_eta.at(type_of_particle)=sparticle->GetEta();
0267           particle_theta.at(type_of_particle)=2*atan(exp(-sparticle->GetEta()));
0268           particle_phi.at(type_of_particle)=sparticle->GetPhi();
0269           particle_pt.at(type_of_particle)=sparticle->GetPt();
0270           particle_p.at(type_of_particle)=sparticle->GetP();
0271           
0272           true_particle_eta.at(type_of_particle)=tparticle->GetEta();
0273           true_particle_theta.at(type_of_particle)=tparticle->GetTheta();
0274           true_particle_phi.at(type_of_particle)=tparticle->GetPhi();
0275           true_particle_p.at(type_of_particle)=tparticle->GetP();
0276           true_particle_pt.at(type_of_particle)=tparticle->GetPt();
0277         }
0278       else if(type_of_particle==4)
0279         {
0280           true_particle_eta.at(type_of_particle)=tparticle->GetEta();
0281           true_particle_theta.at(type_of_particle)=tparticle->GetTheta();
0282           true_particle_phi.at(type_of_particle)=tparticle->GetPhi();
0283           true_particle_p.at(type_of_particle)=tparticle->GetP();
0284           true_particle_pt.at(type_of_particle)=tparticle->GetPt();
0285         }
0286       
0287     } // end loop over particles
0288 
0289       
0290       bool can_we_reconstruct_3_final_state_leptons = true;
0291       if((particle_p.at(0)<1||particle_p.at(1)<1)||particle_p.at(2)<1)
0292     {
0293       can_we_reconstruct_3_final_state_leptons = false;
0294     }
0295       
0296       bool can_we_reconstruct_scattered_lepton = true;
0297       if(particle_p.at(1)<1)
0298     {
0299       can_we_reconstruct_scattered_lepton = false;
0300     }
0301 
0302       bool can_we_reconstruct_decay_particles = true;
0303       if(particle_p.at(0)<1||particle_p.at(2)<1)
0304     {
0305       can_we_reconstruct_decay_particles = false;
0306     }
0307 
0308       // From the eta, phi and p arrays, calculate the invariant mass
0309       // First, try the decay electron and decay positron
0310       if(can_we_reconstruct_3_final_state_leptons)
0311     {
0312       float m = sqrt(2*particle_pt.at(0)*particle_pt.at(2)*(cosh(particle_eta.at(0)-particle_eta.at(2))-cos(particle_phi.at(0)-particle_phi.at(2))));
0313       invariant_reco_decay=m;
0314     }
0315 
0316       // Now try the scattered electron and decay positron
0317       
0318       if(can_we_reconstruct_3_final_state_leptons)
0319     {
0320       float m = sqrt(2*particle_pt.at(1)*particle_pt.at(2)*(cosh(particle_eta.at(1)-particle_eta.at(2))-cos(particle_phi.at(1)-particle_phi.at(2))));
0321       invariant_reco_scatter=m;
0322     }
0323       // cout << particle_theta.at(0) << endl;
0324       // Next, fill in reco J/Psi stuff 
0325       if(can_we_reconstruct_3_final_state_leptons)
0326     {
0327       float sum_px = particle_pt.at(0)*cos(particle_phi.at(0))+particle_pt.at(2)*cos(particle_phi.at(2));
0328       float sum_py = particle_pt.at(0)*sin(particle_phi.at(0))+particle_pt.at(2)*sin(particle_phi.at(2));
0329       float sum_pz = particle_p.at(0)*cos(particle_theta.at(0))+particle_p.at(2)*cos(particle_theta.at(2));
0330       float sum_p = sqrt(sum_px*sum_px+sum_py*sum_py+sum_pz*sum_pz);
0331       jpsi_px_reco=sum_px;
0332       jpsi_py_reco=sum_py;
0333       jpsi_pz_reco=sum_pz;
0334       
0335       
0336       float sum_px_true = true_particle_pt.at(0)*cos(true_particle_phi.at(0))+true_particle_pt.at(2)*cos(true_particle_phi.at(2));
0337       float sum_py_true = true_particle_pt.at(0)*sin(true_particle_phi.at(0))+true_particle_pt.at(2)*sin(true_particle_phi.at(2));
0338       float sum_pz_true = true_particle_p.at(0)*cos(true_particle_theta.at(0))+true_particle_p.at(2)*cos(true_particle_theta.at(2));      
0339       jpsi_px_truth=sum_px_true;
0340       jpsi_py_truth=sum_py_true;
0341       jpsi_pz_truth=sum_pz_true;
0342 
0343       
0344       
0345       // if(sum_p>4.8&&sum_p<5)
0346       //cout << ievent << " - " << true_particle_p.at(4) << " : " << sum_p << " : " << sum_pz << " : " << sum_px << " : " << sum_py << " : " << endl;
0347       //float jpsi_eta = -log(tan( (acos( (particle_p.at(0)*cos(2*atan(exp(-particle_eta.at(0))))+particle_p.at(2)*cos(2*atan(exp(-particle_eta.at(2))))/(sum_p) )))/2));
0348       //float jpsi_eta = -log(tan((acos(particle_p.at(0)*cos(particle_theta.at(0))+particle_p.at(2)*cos(particle_theta.at(2))))
0349       // float jpsi_eta = -log(tan(atan( (particle_pt.at(0)+particle_pt.at(2))/(particle_p.at(0)*cos(particle_theta.at(0))+particle_p.at(2)*cos(particle_theta.at(2))))/2));
0350       //      float jpsi_phi = acos( (particle_pt.at(0)*cos((particle_phi.at(0)))+particle_pt.at(2)*cos((particle_phi.at(2))))/(sqrt(sum_px*sum_px+sum_py*sum_py)) );
0351       float jpsi_eta = -log(tan( (acos(jpsi_pz_reco/sum_p))/2 ) );
0352       float jpsi_phi = atan(jpsi_py_reco/jpsi_px_reco);
0353       if(jpsi_py_reco>0&&jpsi_px_reco<0)
0354         jpsi_phi=3.14159265+jpsi_phi;
0355       else if(jpsi_py_reco<0&&jpsi_px_reco<0)
0356         jpsi_phi=3.14159265+jpsi_phi;
0357       else if(jpsi_py_reco<0&&jpsi_px_reco>0)
0358         jpsi_phi=2*3.14159265+jpsi_phi;
0359         
0360       //cout << true_particle_phi.at(4) << " : " << jpsi_px_truth << " : " << jpsi_py_truth << " : " << atan(jpsi_py_truth/jpsi_px_truth)<< endl;
0361       //cout << true_particle_phi.at(4) << " : " << jpsi_phi<< endl;
0362         //cout << jpsi_phi << " : " << true_particle_phi.at(4) << endl;
0363       particle_p.at(4)=sum_p;
0364       particle_pt.at(4)=sum_p*sin(jpsi_phi);
0365       particle_eta.at(4)=jpsi_eta;
0366       particle_phi.at(4)=jpsi_phi; 
0367       
0368 
0369       // Event Kinematics
0370       
0371       
0372       t = 2*275*true_particle_p.at(3)*(1-cos(2*atan(exp(-true_particle_eta.at(3)))));
0373     }
0374       
0375       de_phi_truth=true_particle_phi.at(0);
0376       de_eta_truth=true_particle_eta.at(0);
0377       de_pt_truth=true_particle_pt.at(0);
0378       de_p_truth=true_particle_p.at(0);
0379       dp_phi_truth=true_particle_phi.at(2);
0380       dp_eta_truth=true_particle_eta.at(2);
0381       dp_pt_truth=true_particle_pt.at(2);
0382       dp_p_truth=true_particle_p.at(2);
0383       sp_phi_truth=true_particle_phi.at(3);
0384       sp_eta_truth=true_particle_eta.at(3);
0385       sp_pt_truth=true_particle_pt.at(3);
0386       sp_p_truth=true_particle_p.at(3);
0387       se_phi_truth=true_particle_phi.at(1);
0388       se_eta_truth=true_particle_eta.at(1);
0389       se_pt_truth=true_particle_pt.at(1);
0390       se_p_truth=true_particle_p.at(1);
0391       jpsi_phi_truth=true_particle_phi.at(4);
0392       jpsi_eta_truth=true_particle_eta.at(4);
0393       jpsi_pt_truth=true_particle_pt.at(4);
0394       jpsi_p_truth=true_particle_p.at(4);
0395       de_phi_reco=particle_phi.at(0);
0396       de_eta_reco=particle_eta.at(0);
0397       de_pt_reco=particle_pt.at(0);
0398       de_p_reco=particle_p.at(0);
0399       dp_phi_reco=particle_phi.at(2);
0400       dp_eta_reco=particle_eta.at(2);
0401       dp_pt_reco=particle_pt.at(2);
0402       dp_p_reco=particle_p.at(2);
0403       sp_phi_reco=particle_phi.at(3);
0404       sp_eta_reco=particle_eta.at(3);
0405       sp_pt_reco=particle_pt.at(3);
0406       sp_p_reco=particle_p.at(3);
0407       se_phi_reco=particle_phi.at(1);
0408       se_eta_reco=particle_eta.at(1);
0409       se_pt_reco=particle_pt.at(1);
0410       se_p_reco=particle_p.at(1);
0411       jpsi_phi_reco=particle_phi.at(4);
0412       jpsi_eta_reco=particle_eta.at(4);
0413       jpsi_pt_reco=particle_pt.at(4);
0414       jpsi_p_reco=particle_p.at(4);
0415       // if(ievent==770)
0416       //    cout << jpsi_p_truth << " : " << jpsi_p_reco << endl;
0417       theTree->Fill();
0418     } // end loop over events
0419   
0420   /////
0421 
0422  
0423 
0424   myfile->cd();
0425   theTree->Write("Tree");
0426   myfile->Write();
0427   myfile->Close();
0428   return 0;
0429 }
0430 
0431 float eta_to_theta(float eta)
0432 {
0433   return 2*atan(exp(-1*eta));
0434 }
0435 
0436 float theta_to_eta(float theta)
0437 {
0438   return -log(tan(theta/2));
0439 }