Upload New File

analysesimuGen.C

+#include "analysesimuGen.h"
+#include "TFile.h"
+#include "TTree.h"
+#include "TChain.h"
+#include "TGraphErrors.h"
+#include "TH1F.h"
+#include "TF1.h"
+#include "TF2.h"
+#include "TTreeReader.h"
+#include "TTreeReaderValue.h"
+#include "TTreeReaderArray.h"
+#include "TH2D.h"
+#include "TLorentzVector.h"
+#include "TVector3.h"
+#include "TMath.h"
+#include "TCanvas.h"
+#include "TH3F.h"
+#include "TString.h"
+#include "reader.h"
+#include "TStyle.h"
+#include "MParticle.h"
+#include "TParticle.h"
+#include "dirent.h"
+#include <iostream>
+
+
+
+using namespace std;
+
+analysesimuGen::analysesimuGen(){
+
+    c1 = new TCanvas("","",800,600);
+    c1->Divide(2,1);
+
+
+    c2 = new TCanvas("Graph errors","Graph errors", 800,600);    
+    c2->Divide(2,2);
+
+     c3 = new TCanvas("cos","cos", 200,10,700,500);
+	 c3->Divide(3,2);
+
+
+    c4 = new TCanvas("Kinematic observables","Kinematic observables",800,600);
+    c4->Divide(2,3);
+
+     c5 = new TCanvas("histo","histo",800,600);
+       c5->Divide(1,2);
+
+    c7 = new TCanvas("test","test",800,600);
+    c7-> Divide(2,1);
+
+    //c6 = new TCanvas("histo","histo",800,600);
+    //c6->Divide(2,2);
+    
+    // bin limits
+
+    binz = new double[nbbinz+1]{0.31,0.38,0.5,0.62,0.8};
+    binx = new double[nbbinx+1]{0.1,0.2,0.3,0.38,1};
+    binp = new double[nbbinp+1]{0.05,0.3,0.45,0.62,1.2};
+    binq = new double[nbbinq+1]{1,2,2.5,3.3,15};
+
+    int a=0;
+
+    for (int i=0; i<nbbinx; i++) {
+        for (int j=0; j<nbbinz; j++){
+            for (int k=0; k<nbbinp; k++){
+                for (int l=0; l<nbbinq; l++){
+
+                    xx[i]=0;
+                    ii[i]=0;
+                    meanbinx[i]=0;
+
+                    cosphi[i][j][k][l]  = new TH1F(Form("h01_%d%d%d%d", i,j,k,l),"",100,-1,1);
+                    cos2phi[i][j][k][l] = new TH1F(Form("h02_%d%d%d%d", i,j,k,l),"",100,-1,1);
+                    cossphi[i][j][k][l] = new TH1F(Form("h0s_%d%d%d%d", i,j,k,l),"",100,-1,1);
+
+                    if (i==0){
+                        graphx[j][k][l] = new TGraphErrors(nbbinx);
+                        graphx[j][k][l]->GetXaxis()->SetTitle("x_B");
+                        graphx[j][k][l]->GetYaxis()->SetTitle("<cos(phi)>");
+                        graphx[j][k][l]->SetTitle(Form("graphx_%d%d%d",j,k,l));
+
+                        graphx2[j][k][l] = new TGraphErrors(nbbinx);
+                        graphx2[j][k][l]->GetXaxis()->SetTitle("x_B");
+                        graphx2[j][k][l]->GetYaxis()->SetTitle("<cos(2phi)>");
+                        graphx2[j][k][l]->SetTitle(Form("graphx2_%d%d%d",j,k,l));
+			
+                    }
+
+                    if (j==0){
+                        graphz[i][k][l] = new TGraphErrors(nbbinz);
+                        graphz[i][k][l]->GetXaxis()->SetTitle("z");
+                        graphz[i][k][l]->GetYaxis()->SetTitle("<cos(phi)>");
+                        graphz[i][k][l]->SetTitle(Form("graphz_%d%d%d",i,k,l));
+
+
+                        graphz2[i][k][l] = new TGraphErrors(nbbinz);
+                        graphz2[i][k][l]->GetXaxis()->SetTitle("z");
+                        graphz2[i][k][l]->GetYaxis()->SetTitle("<cos(2phi)>");
+                        graphz2[i][k][l]->SetTitle(Form("graphz2_%d%d%d",i,k,l));
+			
+                    }
+
+                    if (k==0){
+                        graphp[i][j][l] = new TGraphErrors(nbbinp);
+                        graphp[i][j][l]->GetXaxis()->SetTitle("P_t");
+                        graphp[i][j][l]->GetYaxis()->SetTitle("<cos(phi)>");
+                        graphp[i][j][l]->SetTitle(Form("graphp_%d%d%d",i,j,l));
+
+                        graphp2[i][j][l] = new TGraphErrors(nbbinp);
+                        graphp2[i][j][l]->GetXaxis()->SetTitle("P_t");
+                        graphp2[i][j][l]->GetYaxis()->SetTitle("<cos(2phi)>");
+                        graphp2[i][j][l]->SetTitle(Form("graphp2_%d%d%d",i,j,l));
+			
+			
+                    }
+
+                    if (l==0){
+                        graphq[i][j][k] = new TGraphErrors(nbbinq);
+                        graphq[i][j][k]->GetXaxis()->SetTitle("Q2");
+                        graphq[i][j][k]->GetYaxis()->SetTitle("<cos(phi)>");
+                        graphq[i][j][k]->SetTitle(Form("graphq_%d%d%d",i,j,k));
+
+
+                        graphq2[i][j][k] = new TGraphErrors(nbbinq);
+                        graphq2[i][j][k]->GetXaxis()->SetTitle("Q2");
+                        graphq2[i][j][k]->GetYaxis()->SetTitle("<cos(2phi)>");
+                        graphq2[i][j][k]->SetTitle(Form("graphq2_%d%d%d",i,j,k));
+			
+                    }
+
+                }
+            }
+        }
+    }
+
+
+
+        hQ2 = new TH1F("Q2", "", 100,1,10);
+	hW = new TH1F("W", "", 100,1.8,4);
+	hEcalotot = new TH1F("Ecalotot","Ecalotot",100,0,10);
+	hnpheche = new TH1F("nphe in the cherenkov","nphe in the cherenkov", 100,0,50);
+	hsf = new TH1F("sf","sf",100,0,10);
+	hpz = new TH1F("pz","pz",100,-1.1,1.1);
+	hP = new TH1F("P","P",100,-1.1,1.1);
+	hpid = new TH1F("pid","pid",100,-300,300);
+	hnphe = new TH1F("nphe","nphe",100,-150,200);
+	hEcal = new TH1F("Ecal","Ecal",100,0,800);
+	hducalo = new TH1F("ducalo","ducalo",100,0,10);
+	hy = new TH1F("y","",100,0,0.85);
+	hxB = new TH1F("xB","",100,0,1);
+	hz = new TH1F("z","",100,0,5); 
+	hPt = new TH1F("Pt","",100,0,10);
+	hphi = new TH1F("phi","",100,0,360);
+	hcos = new TH1F("cos phi","cos phi",100,0,1);
+	hcos2 = new TH1F("cos (2phi)","",100,-1,1);
+	hzbin = new TH1F("z","",4,binz);
+	hxbin = new TH1F("x","",4,binx);
+	hQbin = new TH1F("Q2","",4,binq);
+	hPbin = new TH1F("Pt","",4,binp);
+	hMissingM = new TH1F("Missing mass","",100,0,5);
+	hMissingM2 = new TH1F("Missing mass","",100,0.5,2);
+	
+	hnumz = new TH1F("hnumz","",4,binz);
+	hnumQ = new TH1F("hnumQ","",4,binq);
+	hnumP = new TH1F("hnumP","",4,binp);
+	hnumx = new TH1F("hnumx","",4,binx);
+	
+	hnpheVSdet = new TH2F("hnpheVSdet", "hnpheVSdet",100,0,20,100,-150,150);
+	hEcalototVSP = new TH2F("EcalototVSP","EcalototVSP", 100,0,20,100,0,20);
+	hQ2vsW = new TH2F("Q2 vs W", "",100,1,10,100,1.8,4);
+	hQ2vsxB = new TH2F("Q2 vs xB", "",100,1,10,100,0,1);
+	hcosvsz = new TH2F("","",4,0,0.5,4,0,2);
+	hmmvsQ2 = new TH2F("hmmvsW", "hmmvsW",100,0,10,100,0,5);
+    
+Q2=0;
+W =0;
+mu =0;
+y =0;
+xB =0;
+z =0;
+Pt =0;
+l;
+phi =0;
+sign =0;
+cosh=0;
+nbevent=0;
+
+
+Ebeam=10.6;
+M=0.938;
+Q2cut = 1.5;
+Wcut = 2;
+ycut = 0.8;
+zmincut=0.3;
+zmaxcut=0.8;
+Ptcut=1.2;
+Missingm=0;
+ 
+entries=0;
+mean=0;
+error=0;
+
+entries2=0;
+mean2=0;
+error2=0;     
+
+sentries=0;
+smean=0;
+serror=0;
+sigma=0;
+sigma2=0;
+gety=0;
+nbelectron=0;
+nbpion=0;
+nbelectron2=0;
+nbpion2=0;
+nbz=0;
+ nbz2=0;
+ nbmean=0; 
+// part=new MParticle[np];
+ nbread=0;
+
+}
+
+
+analysesimuGen::~analysesimuGen(){
+
+  delete c1,c2,c3,c4,c5,c6;
+    delete binx, binz, binq, binp;
+
+    delete hQ2,hW,hEcalotot,hnpheche,hsf,hpz,hP,hpid,hnphe,hEcal,hducalo,hy,hxB,hz,hPt,hphi,hcos,hcos2,hzbin,hxbin,hQbin,hPbin,hnumz,hnumQ,hnumP,hnumx,hnpheVSdet,hEcalototVSP,hQ2vsW,hQ2vsxB,hcosvsz;
+
+    
+    //   delete part;
+
+    for (int i=0; i<nbbinx; i++) {
+       for (int j=0; j<nbbinz; j++){
+           for (int k=0; k<nbbinp; k++){
+               for (int l=0; l<nbbinq; l++){
+
+		 delete cosphi[i][j][k][l],cos2phi[i][j][k][l],cossphi[i][j][k][l],graphx[j][k][l],graphz[i][k][l], graphp[i][j][l], graphq[i][j][k],graphx2[j][k][l],graphz2[i][k][l], graphp2[i][j][l], graphq2[i][j][k];
+
+		   }
+	   }
+       }
+    }
+
+}
+
+void analysesimuGen::open_file(TString nameFiles){
+
+  //	nameFiles="/home/caudron/Documents/Analyse/Filtered.hipo";
+
+
+
+	//hipo reader
+	//Define all the stuff you need to read a file
+
+	
+       
+	reader.open(nameFiles);
+
+	reader.readDictionary(factory);
+	factory.show();
+
+	
+	// create the bank container
+	// you can add new bank you want to read
+
+
+	// Create canvas and histograms
+
+
+	float nbevent=0;
+
+
+
+}
+
+int analysesimuGen::read_eventdata(){
+  //  cout<<"read"<<endl;
+
+  //  cout<<nbread<<endl;
+  // cout<<reader.next()<<endl;
+   if (reader.next()==0){
+     return 1;}
+
+	
+   /*	hipo::bank EVENT(factory.getSchema("REC::Event"));
+	hipo::bank PART(factory.getSchema("REC::Particle"));
+	hipo::bank SCIN(factory.getSchema("REC::Scintillator"));
+	hipo::bank CHE(factory.getSchema("REC::Cherenkov"));
+	hipo::bank CALO(factory.getSchema("REC::Calorimeter"));
+	hipo::bank RUN(factory.getSchema("RUN::config"));*/
+	hipo::bank MCPART(factory.getSchema("MC::Particle"));
+	hipo::bank MCEVENT(factory.getSchema("MC::Event"));
+	
+Q2=0;
+W =0;
+mu =0;
+y =0;
+xB =0;
+z =0;
+Pt =0;
+l;
+phi =0;
+sign =0;
+cosh=0;
+  
+		//read the event
+		reader.read(event);
+
+		//get the bank in the container
+		//this has to be done to read the data
+		event.getStructure(MCPART);
+		event.getStructure(MCEVENT);
+		/*	event.getStructure(PART);
+		event.getStructure(SCIN);
+		event.getStructure(CHE);
+		event.getStructure(CALO);
+		event.getStructure(EVENT);*/
+
+		//show Particle bank use to debug
+		//PART.show();
+
+
+		float Ebeam=10.6;
+		float M=0.938;
+		float Q2cut = 1.5;
+		float Wcut = 2;
+		float ycut = 0.8;
+		float zmincut=0.3;
+		float zmaxcut=0.8;
+		float Ptcut=1.2;
+
+		MParticle Electron;
+		MParticle proton;
+		MParticle beam;
+		MParticle electron;
+		MParticle pionm;
+		MParticle pionp;
+		MParticle photon;
+		MParticle MissingM;
+
+
+		proton.Vector.SetPxPyPzE(0,0,0,M);
+		beam.Vector.SetPxPyPzE(0,0,Ebeam,Ebeam);
+		MissingM.Vector.SetPxPyPzE(0,0,0,Missingm);
+
+
+
+		//get the number of particle
+
+		
+	        np=MCPART.getRows();
+	       	class MParticule;
+
+
+		int nevent = MCEVENT.getRows();
+		/*	int nscin = SCIN.getRows();
+		int nche=CHE.getRows();
+		int ncalo=CALO.getRows();*/
+
+
+		MParticle *part=new MParticle[np];
+
+
+		
+		//loop through particles
+		for( int i=0; i < np; i++ ){
+		  
+		  //  cout<<"i       "<<i<<endl;
+
+			//get x,y,z,pid component of the i-th particle from the Particle bank
+			//get nb photoe- in Cherenkov bank
+			//get energy in  calorimeter bank
+
+			float px = MCPART.getFloat("px",i);
+			float py = MCPART.getFloat("py",i);
+			float pz = MCPART.getFloat("pz",i);
+			float vx = MCPART.getFloat("vx",i);
+			float vy = MCPART.getFloat("vy",i);
+			float vz = MCPART.getFloat("vz",i);
+			int pid = MCPART.getInt("pid",i);
+			//	float status=MCPART.getFloat("status",i);
+			//	float chi2pid=MCPART.getFloat("chi2pid",i);
+
+			part[i]=MParticle(pid,0,0,0,0,0,px,py,pz,0,vx,vy,vz,0);
+
+			//	cout << pid<<endl;
+
+
+
+
+			float P = pow(part[i].px_part,2)+pow(part[i].py_part,2)+pow(part[i].pz_part,2);
+			float Efine = sqrt(P);
+			float Efinp= sqrt(P+pow((0.139),2));
+        
+
+
+			if (part[i].GetPdgCode()==11) {
+
+				Electron.Vector.SetPxPyPzE(part[i].px_part,part[i].py_part,part[i].pz_part,Efine);
+
+				nbelectron+=1;
+
+				nbelectron2++;
+
+				//	cout << "e        "<< Efin<<endl;
+
+				
+				//    	cout<<"nbelectron      "<<nbelectron<<endl;
+
+				
+				photon.Vector.SetPxPyPzE(((beam.Vector).Px()-(Electron.Vector).Px()),((beam.Vector).Py()-(Electron.Vector).Py()),((beam.Vector).Pz()-(Electron.Vector).Pz()),((beam.Vector).E()-(Electron.Vector).E()));
+			}
+
+
+			if (part[i].GetPdgCode()==211){
+
+			  	pionp.Vector.SetPxPyPzE(part[i].px_part,part[i].py_part,part[i].pz_part,Efinp);
+
+			  
+				//	pionp.Vector.SetPxPyPzE(part[i].px_part,part[i].py_part,part[i].pz_part,0.139);
+				
+				nbpion+=1;
+				
+
+				nbpion2++;
+				
+				//	cout<<"nbpion        "<<nbpion<<endl;
+
+				//	cout << "p      "<<Efinp<<endl;
+				
+			}
+
+			
+
+			//	cout << "e-  " << nbelectron <<"       pion         " << nbpion << endl;
+
+			
+		}
+
+		//	cout<< "e-       "<< nbelectron2<<"      pion     "<<nbpion2<<endl;
+		//cout << "e-  " << nbelectron <<"       pion         " << nbpion << endl;
+		
+		if (nbelectron>0 && nbpion>0){
+		  nbread++;
+		  nbelectron=0;
+		  nbpion=0;
+		  //	  cout<<nbread<<endl;
+			  //	cout<< "e-2      "<< nbelectron2<<"       pion2     "<<nbpion2<<endl;
+
+		  //	  cout<<Q2<<endl;
+Q2 = (-1*(beam.Vector-Electron.Vector).M2());
+W = sqrt(2*M*((beam.Vector-Electron.Vector).E())+pow(M,2)-Q2);
+mu = (beam.Vector-Electron.Vector).E();
+y =(mu)/((beam.Vector).E());
+xB =(Q2/(2*((proton.Vector).Dot(beam.Vector-Electron.Vector))));
+z =((pionp.Vector).E())/mu;
+
+ nbz++;
+ /* cout <<"z          "<<z<<endl;
+    cout <<"nbz        "<<nbz<<endl;*/
+ 
+Pt = (pionp.Vector).Perp((Electron.Vector).Vect()-(beam.Vector.Vect()));
+l=((Electron.Vector).Vect()).Cross((beam.Vector).Vect());
+phi = (l.Angle(((photon.Vector).Vect()).Cross((pionp.Vector).Vect())))*180/TMath::Pi();
+sign =(l.Cross(((photon.Vector).Vect()).Cross((pionp.Vector).Vect()))).Z();
+
+
+Missingm=((beam.Vector)+(proton.Vector)-(electron.Vector)-(pionp.Vector)).M();
+
+
+	}
+ 
+/*if (part[i].GetPdgCode()==211){
+  //cout<<Q2<<W<<mu<<y<<xB<<z<<Pt<<phi<<endl;
+
+
+  }*/
+			if (sign<0){ phi=360-phi; }
+
+	if (Q2>Q2cut && W>Wcut && y<ycut && z<zmaxcut && z>zmincut && Pt<Ptcut ){
+			for (int i=0;i<nbbinx;i++){
+				if(xB>binx[i] && xB<binx[i+1]){
+					for (int j=0;j<nbbinz;j++){
+						if(z>binz[j] && z<binz[j+1]){
+							for (int k=0;k<nbbinp;k++){
+								if(Pt>binp[k] && Pt<binp[k+1]){
+									for (int l=0;l<nbbinq;l++){
+										if(Q2>binq[l] && Q2<binq[l+1]){
+
+										  //    xx[i]+=xB;
+										  //	    ii[i]++;
+
+										  
+										    //     cout << "i :    " << ii[i] << "   xx :   " << xx[i]<< "   xB :   "<< xB <<endl; 
+											cos2phi[i][j][k][l]->Fill(cos(2*phi));
+											cosphi[i][j][k][l]->Fill(cos(phi));
+
+											//	cout << cos(2*phi)<<endl;
+											cossphi[i][j][k][l]->Fill(pow(cos(2*phi),2));
+
+										}
+									}
+								}
+							}
+						}
+					}
+				}
+			}
+
+	}
+
+
+		
+				nbpion=0;
+			nbelectron=0;	
+
+		return 0;
+	}
+
+/*void analyse::fill_histo(){
+
+
+  		for( int i=0; i < np; i++ ){
+
+		  //	  cout<<Q2<<endl;
+
+
+
+			
+		  }*/
+
+
+
+void analysesimuGen::fill_graph(){
+
+  //	for( int i=0; i < np; i++ ){
+		
+	entries=0;
+	mean=0;
+	error=0;
+	
+	entries2=0;
+	mean2=0;
+	error2=0;     
+
+	sentries=0;
+	smean=0;
+	serror=0;
+	sigma=0;
+	sigma2=0;
+	gety=0;
+
+
+
+			for (int i=0;i<nbbinx;i++){
+				if(xB>binx[i] && xB<binx[i+1]){
+					for (int j=0;j<nbbinz;j++){
+						if(z>binz[j] && z<binz[j+1]){
+							for (int k=0;k<nbbinp;k++){
+								if(Pt>binp[k] && Pt<binp[k+1]){
+									for (int l=0;l<nbbinq;l++){
+										if(Q2>binq[l] && Q2<binq[l+1]){
+
+
+
+
+	
+
+										  /*   for (int i=0;i<nbbinx;i++){
+        for (int j=0;j<nbbinz;j++){
+            for (int k=0;k<nbbinp;k++){
+	    for (int l=0;l<nbbinq;l++){*/
+
+
+		  Int_t axis;
+
+		  
+                    mean=cosphi[i][j][k][l]->GetMean(axis=1);
+
+		      nbmean++;
+		    
+		    //    cout <<mean<<endl;
+
+
+		    
+                    error=cosphi[i][j][k][l]->GetMeanError();
+                    entries=cosphi[i][j][k][l]->GetEntries();
+
+                    mean2=cos2phi[i][j][k][l]->GetMean();
+                    error2=cos2phi[i][j][k][l]->GetMeanError();
+                    entries2=cos2phi[i][j][k][l]->GetEntries();
+
+                    /*	//	error=mean/sqrt(entries);
+                    //	gety = cosphi[i][j][k][l]->GetBinContent();
+
+                    //			sigma = pow((gety-mean),2)/(entries-1);
+                    //	sigma2 = pow((gety-mean2),2)/(entries2-1);
+
+                    //	cout << "cos 2phi =    " << cos(2*phi) << "     mean =    "<< mean << "      entries =     " << entries << endl;
+
+                    smean=cossphi[i][j][k][l]->GetMean();
+                    sentries=cossphi[i][j][k][l]->GetEntries();
+                    //serror=smean/sqrt(sentries);
+                    serror = sqrt(2*smean/sentries);*/
+
+                    graphx[j][k][l]->SetPoint(i, binx[i]+(binx[i+1]-binx[i])/2, mean);
+                    graphx[j][k][l]->SetPointError(i, 0., error);
+
+                    graphx2[j][k][l]->SetPoint(i, binx[i]+(binx