[Hps-analysis] mad5VSdata!
Luca Colaneri
luca.colaneri at roma2.infn.it
Tue Jul 26 09:54:56 EDT 2016
Good news everyone!
mad5 recon data looks good! pdf in attachment
I also attach the piece of code I've been using so far to read recon
data (both exp and MC) using the slcio c++ api.
Bests
L.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: mad5VSdata.pdf
Type: application/pdf
Size: 27228 bytes
Desc: not available
URL: <https://mailman.jlab.org/pipermail/hps-analysis/attachments/20160726/5cacf922/attachment-0001.pdf>
-------------- next part --------------
#include <TH1D.h>
#include <TH2D.h>
#include <TFile.h>
#include <TMath.h>
#include <TCanvas.h>
#include <TChain.h>
#include <TObjArray.h>
//================ LCIO ==================
#include "IO/LCReader.h"
#include "IOIMPL/LCFactory.h"
#include <IMPL/LCCollectionVec.h>
#include "EVENT/LCCollection.h"
//#include "EVENT/Cluster.h"
//#include <EVENT/CalorimeterHit.h>
#include <EVENT/ReconstructedParticle.h>
//#include <ECalCluster.h>
#include <cmath>
#include <fstream>
#include <iostream>
#include <sstream>
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <fcntl.h>
#include <unistd.h>
#include <string>
#include <TApplication.h>
using namespace std;
int main( int argc, const char *argv[] )
{
double pi=acos(-1);
cout<<"So it begins.."<<endl;
IO::LCReader* lcReader = IOIMPL::LCFactory::getInstance()->createLCReader();
long long ind_global = 0;
TH1D *inva=new TH1D("inva","inva",70,0,140);
///istogrammi stati finali
TH1D *enem=new TH1D("Em","Em",100,0,1.05);
TH1D *enep=new TH1D("Ep","Ep",100,0,1.05);
TH2D *ene2d=new TH2D("Ep VS Em","Ep VS Em",100,0,1,100,0,1);
TH1D *eneSum=new TH1D("Ep+Em","Ep+Em",100,0.1,1.2);
//istogrammi angoli
TH1D *thetap= new TH1D("thetap","thetap",100,0,0.4);
TH1D *thetam= new TH1D("thetam","thetam",100,0,0.4);
//correlazioni v2angoli
TH2D *thpV2=new TH2D("thpV2","thpV2",100,0,0.3,140,0,140);
TH2D *thmV2=new TH2D("v2 VS teta","v2 VS teta",100,0,0.3,140,0,140);
TH1D *thp=new TH1D("theta pos.", "theta pos.", 100, 0.01, 0.2);
TH1D *thm=new TH1D("theta el.", "theta el.", 100, 0.01,0.2);
//angolo interno
TH1D *intang=new TH1D("intang","intanf",100,0,0.3);
TH2D *eneang=new TH2D("eneang","eneang",100,0.1,1.2,100,0,0.3);
/////ANGOLI BERA
TH2D *beralpM=new TH2D("Positron Acceptance","Positron Acceptance",100,pi/2-0.15,pi/2 + 0.15,100,-0.1,0.1);
TH2D *beralmM=new TH2D("Electron1 acceptance","Electron1 Acceptance",100,pi/2-0.15,pi/2 + 0.15,100,-0.1,0.1);
//vertici Positrone
TH1D *xp=new TH1D("x pos","x pos",100, -10,10);
TH1D *yp=new TH1D("y pos","y pos",100, -10,10);
TH1D *zp=new TH1D("z pos","z pos",100, -10,10);
///ene pos up and down
///ene pos up and down
TH1D *eneposup=new TH1D("Ep up","Ep up",100,0,1.05);
TH1D *eneposdown=new TH1D("Ep down","Ep down",100,0,1.05);
//theta pos up and down
TH1D *thetaposup=new TH1D("thetap up","theta up",100,0,0.2);
TH1D *thetaposdown=new TH1D("thetap down","thetam down",100,0,0.2);
///ene pos up and down
TH1D *eneeleup=new TH1D("Em up","Em up",100,0,1.05);
TH1D *eneeledown=new TH1D("Em down","Em down",100,0,1.05);
//theta pos up and down
TH1D *thetamup=new TH1D("thetam up","thetam up",100,0,0.2);
TH1D *thetamdown=new TH1D("thetam down","thetam down",100,0,0.2);
EVENT::LCEvent* ev = NULL;
TFile *file_out = new TFile("invaRUN.root", "Recreate");
//string file=argv[1];
string RecPartCol="TargetConstrainedV0Candidates";
double Eplus;
double Eminus;
const double *pp;
const double *pm;
double maxm=0;
double maxp=0;
double minm=100;
double minp=100;
const float *vert;
//std::string word = argv[1];
// ============== open input files ==================
for(int f=1;f<2;f++){
std::ostringstream conv_setting;
conv_setting << f;
string isetting=conv_setting.str(); //convert the for index into a string
string file="hps_005772.0_recon_R3.8.slcio";
lcReader->open(file); ///open file
// int nev = lcReader->getNumberOfEvents();
//cout<<"file "<<argv<<endl;
int ev_number = 0;
while( (ev = lcReader->readNextEvent()) != 0)
{
if( ind_global%100 == 0 )
{
cout.flush()<<"Processed "<<ind_global<<"\r";
}
IMPL::LCCollectionVec *particles = (IMPL::LCCollectionVec*)ev->getCollection(RecPartCol);
int n_particles = particles->getNumberOfElements();
for( int icl = 0; icl < n_particles; icl++ )
{
EVENT::ReconstructedParticle *particle = (EVENT::ReconstructedParticle*)particles->getElementAt(icl);
EVENT::ReconstructedParticleVec candi = (EVENT::ReconstructedParticleVec)particle->getParticles();
EVENT::ReconstructedParticle *positron= NULL;
EVENT::ReconstructedParticle *electron= NULL;
int numpart= candi.size();
int track_type_e;
for(int i=0; i<numpart;i++){
EVENT::ReconstructedParticle *party = (EVENT::ReconstructedParticle*)candi.at(i);
if(party->getCharge()>0 ){positron=party; }
if(party->getCharge()<0 ){electron=party; }
}
Eplus=positron->getEnergy();
Eminus=electron->getEnergy();
pp=positron->getMomentum();
pm=electron->getMomentum();
//angolibera
double lmv=sqrt(pm[0]*pm[0]+pm[2]*pm[2]+pm[1]*pm[1]);
double thblm=asin(pm[1]/lmv);
double alphalm=atan2(pm[2],pm[0]);
double thblmC=sqrt(thblm*thblm);
double lmp=sqrt(pp[1]*pp[1]+pp[2]*pp[2]+pp[0]*pp[0]);
double thblp=asin(pp[1]/lmp);
double alphalp=atan2(pp[2],pp[0]);
double thblpC=sqrt(thblp*thblp);
double thpp=acos(pp[2]/lmp);
double thmm=acos(pm[2]/lmv);
double thetaint;
thetaint=acos((pm[0]*pp[0]+pm[1]*pp[1]+pm[2]*pp[2])/(lmp*lmv));
//thblmC >=0.015 && thblmC <= 0.09 && thblpC >= 0.015 && thblpC <= 0.09 &&
/// IN THIS IF STATEMENT ALL THE CUTS ARE DEFINED
if( thblmC >=0.015 && thblpC >= 0.015 && thblmC <= 0.055 && thblpC <= 0.055
// ( (thblm >= 0.02 && thblm <= 0.055) || (thblm >=-0.055 && thblm <= -0.02) ) && ( (thblp >= 0.02 && thblp <= 0.055) || (thblp >=-0.055 && thblp <= -0.02) )
&& pp[1]*pm[1] <0 && Eplus>0.05 /*&& Eplus<0.9 && Eminus <0.9*/
&& Eminus>0.05 && (Eplus+Eminus)>=0.2 && (Eminus + Eplus)<1.05
&& positron->getTracks().at(0)->getType()>=32 && electron->getTracks().at(0)->getType()>=32 /*&& thetaint<0.05*/){
if(pp[1]>0){
eneposup->Fill(Eplus);
thetaposup->Fill(acos(pp[2]/(sqrt(pp[0]*pp[0]+pp[1]*pp[1]+pp[2]*pp[2]) )));
}
else{
eneposdown->Fill(Eplus);
thetaposdown->Fill( acos(pp[2]/(sqrt(pp[0]*pp[0]+pp[1]*pp[1]+pp[2]*pp[2]) )) );
}
if(pm[1]>0){
eneeleup->Fill(Eminus);
thetamup->Fill(acos(pm[2]/(sqrt(pm[0]*pm[0]+pm[1]*pm[1]+pm[2]*pm[2]) )));
}
else{
eneeledown->Fill(Eminus);
thetamdown->Fill(acos(pm[2]/(sqrt(pm[0]*pm[0]+pm[1]*pm[1]+pm[2]*pm[2]) )));
}
if(thblmC>maxm){maxm=thblmC;}
if(thblmC<minm){minm=thblmC;}
if(thblpC>maxp){maxp=thblpC;}
if(thblpC<minp){minp=thblpC;}
thp->Fill(thpp);
thm->Fill(thmm);
inva->Fill(particle->getMass()*1000 );
enep->Fill(Eplus);
enem->Fill(Eminus);
ene2d->Fill(Eminus,Eplus);
thetap->Fill( acos(pp[2]/(sqrt(pp[0]*pp[0]+pp[1]*pp[1]+pp[2]*pp[2]) )));
thetam->Fill( acos(pm[2]/(sqrt(pm[0]*pm[0]+pm[1]*pm[1]+pm[2]*pm[2]) )));
thpV2->Fill(acos(pp[2]/(sqrt(pp[0]*pp[0]+pp[1]*pp[1]+pp[2]*pp[2]) )), particle->getMass()*1000 );
thmV2->Fill(acos(pm[2]/(sqrt(pm[0]*pm[0]+pm[1]*pm[1]+pm[2]*pm[2]) )), particle->getMass()*1000 );
beralpM->Fill(alphalp,thblp);
beralmM->Fill(alphalm,thblm);
eneSum->Fill(Eplus+Eminus);
intang->Fill(thetaint);
eneang->Fill(Eplus+Eminus,thetaint);
}
}
}
}//fine for
eneposup->Write();
eneposdown->Write();
thetaposup->Write();
thetaposdown->Write();
eneeleup->Write();
eneeledown->Write();
thetamup->Write();
thetamdown->Write();
TCanvas *thcan=new TCanvas ("dd","dd");
thcan->Divide(1,2);
thcan->cd(1);
thp->GetXaxis()->SetTitle("theta positron");
thp->Draw();
thp->Write();
thcan->cd(2);
thm->GetXaxis()->SetTitle("theta electron");
thm->Draw();
thm->Write();
thcan->Print("inva.pdf(");
TCanvas *pospcan=new TCanvas("pos","pos");
pospcan->Divide(2,2);
pospcan->cd(1);
xp->Draw();
xp->Write();
pospcan->cd(2);
yp->Draw();
yp->Write();
pospcan->cd(3);
zp->Draw();
zp->Write();
pospcan->Print("inva.pdf");
TCanvas *mycan=new TCanvas("c","c");
inva->GetXaxis()->SetTitle("Inv.Mass (MeV)");
inva->Draw();
inva->Write();
mycan->Print("inva.pdf");
TCanvas *enecan=new TCanvas("c2","c2");
enecan->Divide(2,2);
enecan->cd(1);
enep->GetXaxis()->SetTitle("Ep (GeV)");
enep->Draw();
enep->Write();
enecan->cd(2);
enem->GetXaxis()->SetTitle("Em (GeV)");
enem->Draw();
enem->Write();
enecan->cd(3);
ene2d->GetXaxis()->SetTitle("Em (GeV)");
ene2d->GetYaxis()->SetTitle("Ep (GeV)");
ene2d->Draw("colz");
ene2d->Write();
enecan->cd(4);
eneSum->GetXaxis()->SetTitle("Ep+Em (GeV)");
eneSum->Draw();
eneSum->Write();
enecan->Print("inva.pdf");
TCanvas *angcan= new TCanvas("c3","c3");
angcan->Divide(2,1);
angcan->cd(1)->SetLogy();
thetap->GetXaxis()->SetTitle("theta p (rad)");
thetap->Draw();
thetap->Write();
angcan->cd(2)->SetLogy();
thetam->GetXaxis()->SetTitle("theta m (rad)");
thetam->Draw();
thetam->Write();
angcan->Print("inva.pdf");
TCanvas *cacca=new TCanvas("ca","ca");
thpV2->Draw("colz");
thpV2->Write();
cacca->Print("inva.pdf");
TCanvas *cacca2=new TCanvas("ca2","ca2");
thmV2->Draw("colz");
thmV2->Write();
cacca2->Print("inva.pdf");
TCanvas *beracanM=new TCanvas("beracanM","beracanM");
beracanM->Divide(2,2);
beracanM->cd(1);
beralpM->GetXaxis()->SetTitle("alpha p (rad)");
beralpM->GetYaxis()->SetTitle("beta p (rad) ");
beralpM->Draw("colz");
beralpM->Write();
beracanM->cd(2);
beralmM->GetXaxis()->SetTitle("alpha e1");
beralmM->GetYaxis()->SetTitle("beta e1");
beralmM->Draw("COLZ");
beralmM->Write();
beracanM->Print("inva.pdf");
TCanvas *intangCan=new TCanvas("1","1");
intangCan->Divide(2,1);
intangCan->cd(1);
intang->GetXaxis()->SetTitle("PpPm angle");
intang->Draw();
intang->Write();
intangCan->cd(2);
eneang->GetYaxis()->SetTitle("PpPm angle (rad)");
eneang->GetXaxis()->SetTitle("Esum");
eneang->Draw("COLZ");
eneang->Write();
intangCan->Print("inva.pdf)");
cout<<"massimo elettrone= "<<maxm<<endl;
cout<<"minimo elettrone= "<<minm<<endl;
cout<<"massimo positrone= "<<maxp<<endl;
cout<<"minimo positrone= "<<minp<<endl;
return 0;
}
More information about the Hps-analysis
mailing list