More

More Examples

This example shows some quick steps in an analysis with some basic histograming.

Histograming Examples

const Float_t vlight=29.9792;
const Float_t xmpi=0.139570;
const Float_t xmkp=0.493677;
const Float_t xmk0=0.497614;
const Float_t xmpr=0.938272;
const Float_t xmlam=1.115683;
const Float_t xmsigp=1.18937;
const Float_t xmsig0=1.19264;
const Float_t xmsigm=1.19745;
const Float_t xmksp=0.89166;
const Float_t cutk0=0.015;
const Float_t cutksp=0.08;

float pg_corr(float);
int ibadpd(int,int);
float xmm3(int,int,int,float,float,float,float);
float invm(int,int,int,float,float,float);

TChain *EV=new TChain("h10");
g11 *test = new g11(EV);
extern"C"{
   void fortfunc_(int *ii, float *ff);
}

int main(int argc, char **argv)
{
   char *outFileName = "OutPutFile.root";
   TROOT troot();  
   TFile outFile(outFileName,"recreate");

int ii=5;
   float ff=5.5;
   fortfunc_(&ii, &ff);

TString dname="/w/work5713/tangwei/learn_root/dat1";// Calculate the data 
   cout<<dname<<endl;
   void *dirp = gSystem->OpenDirectory(dname);     
   const char * afile;  
   while((afile = gSystem->GetDirEntry(dirp))){  
      TString fname=afile;
      if(fname.Contains(".root")){ 
         cout<<fname<<endl;
         gSystem->ChangeDirectory(dname);
         EV->Add(fname);
      } 
   }

TH1F *TOF1=new TH1F("TOF1","TOF sector 1",100,0.,50.);  
   TH1F *TOF2=new TH1F("TOF2","TOF sector 2",100,0.,50.);  
   TH1F *TOF3=new TH1F("TOF3","TOF sector 3",100,0.,50.);  
   TH1F *TOF4=new TH1F("TOF4","TOF sector 4",100,0.,50.);  
   TH1F *TOF5=new TH1F("TOF5","TOF sector 5",100,0.,50.);  
   TH1F *TOF6=new TH1F("TOF6","TOF sector 6",100,0.,50.);  
   TH1F *XMK01=new TH1F("XMK01","XMK01",100,0.45,0.55);  
   TH1F *XMK02=new TH1F("XMK02","XMK02",100,0.45,0.55);  
   TH1F *XKSPA=new TH1F("XKSP_A","XKSP_A",100,0.3,1.5);  
   TH1F *XKSPB=new TH1F("XKSP_B","XKSP_B",100,0.3,1.5);  
   TH1F *XKSPC=new TH1F("XKSP_C","XKSP_C",100,0.3,1.5);  
   TH1F *XMMA=new TH1F("XMM_A","XMM_A",100,1.0,1.5);  
   TH1F *XMMB=new TH1F("XMM_B","XMM_B",100,1.0,1.5);  
   TH1F *XMMC=new TH1F("XMM_C","XMM_C",100,1.0,1.5);

TH2F *DTOF1=new TH2F("DTOF1","p vs dtof",100,0,3,100,-10.,10.);  
   TH2F *DTOF2=new TH2F("DTOF2","p vs dtof",100,0,3,100,-10.,10.);

//  TH1F *mppc[20];
   //   for(Int_t i=0;i<20;i++){
   //    TString hn=Form("%d_MPPC_%d",atoi(fn),i+1);
   //    mppc[i]= new TH1F(hn,hn,100,0,3);
   //   }

TObjArray *Hlist = new TObjArray();

TH1F *Lambda_A[9][40];
   for(Int_t i=0;i<40;i++){
      for(Int_t j=0;j<9;j++){
	Lambda_A[j][i]=new TH1F(Form("Lambda_A_%d_%d", j, i),"Lambda_A",100,1.0,1.5);
	Lambda_A[j][i]->SetDirectory(0);
        Hlist->Add(Lambda_A[j][i]);
      }
   }

// TH1F *Lambda_B[9][40];
   //   for(Int_t i=0;i<40;i++){
   //      for(Int_t j=0;j<9;j++){
   //	Lambda_B[j][i]=new TH1F("Lambda_B","Lambda_B",100,1.0,1.5);
   //	Lambda_B[j][i]->SetDirectory(0);
   //      }
   //   }

//   TH1F *Lambda_C[9][40];
   //   for(Int_t i=0;i<40;i++){
   //      for(Int_t j=0;j<9;j++){
   //	Lambda_C[j][i]=new TH1F("Lambda_C","Lambda_C",100,1.0,1.5);
   //	Lambda_C[j][i]->SetDirectory(0);
   //      }
   //   }

Int_t nentries = Int_t(EV->GetEntries());
   for(Int_t iev=0;iev<nentries;iev++){

Int_t  jj,kk;
      Int_t  i1,i2,i3;
      Int_t  ipip[5],ipim[5];
      Int_t  ipip1[5],ipim1[5];
      Int_t  ngam,npip,npim,npip1,npim1;
      Int_t  qq,iselect;
      Int_t  ibad,isector;
      Int_t  xvkey1,xvkey2;
      Int_t  mm,nn;

Float_t pp;
      Float_t egam;
      Float_t betap,vtime,tdif,tdif1;
      Float_t kdif11,kdif12,kdif21,kdif22;
      Float_t dif1,dif2;
      Float_t xt1,xt2,xt3;
      Float_t xp1,xp2,xp3;
      Float_t xz1,xz2,xz3;
      Float_t xtp1,xtp2,xtp3;
      Float_t p_x,p_y,p_z,e,pho4;
      Float_t xmk01,xmk02,xksp,xmm,xang;
      Float_t ecm,beta,gam;
      Float_t pkpzcm,pkpcm;
		 
      TLorentzVector P4pim,P4pip1,P4pip2;
      TLorentzVector P4target,P4pho;
      TLorentzVector P4k01,P4k02;  
      TLorentzVector P4ksp,P4tot;

test->GetEntry(iev);

ngam =0;
      egam =0.;

for(Int_t k=0;k<5;k++){
	 ipip[k]=0;
	 ipim[k]=0;
	 ipip1[k]=0;
	 ipim1[k]=0;
      }

for(Int_t j=0; j<test->taggoodhit; j++){
	 npip=0;
         npim=0;
	 Lambda_A[5][23]->Fill(xmlam);
	 
         for(Int_t i=0;i<test->gpart;i++){
	    if(test->stat[i]>0 && test->sc[i]>0){
               jj=test->sc_sect[test->sc[i]-1];
               kk=test->sc_pd[test->sc[i]-1];
               if(jj == 1){
	          TOF1->Fill((float)kk);
               }     
               else if(jj == 2){
    	          TOF2->Fill((float)kk);         
               }
               else if(jj == 3){
    	          TOF3->Fill((float)kk);         
               }
               else if(jj == 4){
    	          TOF4->Fill((float)kk);         
               }
               else if(jj == 5){
    	          TOF5->Fill((float)kk);         
               }
               else if(jj == 6) {
    	          TOF6->Fill((float)kk);         
               }
	       pp=test->p[i];
	       ibad=ibadpd(jj,kk);
	       if(ibad == 0){
	          betap=pp/sqrt(pp*pp+xmpi*xmpi);
	          vtime=test->vertex_time[j];
	          tdif1=test->sc_t[test->sc[i]-1]-vtime;
	          tdif=tdif1-test->sc_r[test->sc[i]-1]/(vlight*betap);
	          DTOF1->Fill(pp,tdif);
	          if(fabs(tdif) < 1.0){
		     DTOF2->Fill(pp,tdif);
		     if(qq > 0){
		        npip=npip+1;
		        ipip[npip]=i;
		     }
		     if(qq < 0){
		        npim=npim+1;
		        ipim[npim]=i;
		     }
	          }
	       }//ibad
            }//stat and sc
	    if(npip==2 && npim >0 && npim < 3){
	       ngam=ngam+1;
	       egam=test->tag_energy[j];
	       ipip1[0]=ipip[0];
	       ipip1[1]=ipip[1];
	       ipim1[0]=ipim[0];
	       ipim1[1]=ipim[1];
	       npip1=npip;
	       npim1=npim;
	    }
         }//gpart
      }//taggoodhit

iselect=0;
	 
      if(ngam > 0){
	 iselect=1;
	 i1=ipip1[0];
	 i2=ipip1[1];
	 i3=ipim1[0];
	 if(npim1 == 2){
	    kdif11=fabs(xmk0-invm(i1,0,ipim1[0],xmpi,0.,xmpi));
 	    kdif12=fabs(xmk0-invm(i1,0,ipim1[1],xmpi,0.,xmpi));
	    kdif11=fabs(xmk0-invm(0,i2,ipim1[0],0.,xmpi,xmpi));
	    kdif11=fabs(xmk0-invm(0,i2,ipim1[1],0.,xmpi,xmpi));
	    if(kdif11 < kdif21){
	       dif1=kdif11;
	    }	    
	    else{
	       dif1=kdif21;
	    }
	    if(kdif12 < kdif22){
	       dif2=kdif12;
	    }	    
	    else{
	       dif2=kdif22;
	    }	       
	    if(dif1 > dif2){
	       i3=ipim1[1];
	    }
	 }
      }

if(iselect == 1){

xt1=acos(test->cz[i1]);
	xt2=acos(test->cz[i2]);
	xt3=acos(test->cz[i3]);

if(xt1 > 0.135 && xt2 > 0.135 && xt3 > 0.235){
	   
	   xp1=fabs(atan2(test->cz[i1],test->cz[i1]));
	   isector=(xp1+0.52)/1.05;
	   xp1=fabs(xp1-1.05*isector);
	   xtp1=2.55*(1.0-exp(-14.0*xt1))-2.1;
	  
	   xp2=fabs(atan2(test->cz[i2],test->cz[i2]));
	   isector=(xp2+0.52)/1.05;
	   xp2=fabs(xp2-1.05*isector);
	   xtp2=2.55*(1.0-exp(-14.0*xt2))-2.1;

xp3=fabs(atan2(test->cz[i3],test->cz[i3]));
	   isector=(xp3+0.52)/1.05;
	   xp3=fabs(xp3-1.05*isector);
	   xtp3=2.55*(1.0-exp(-14.0*(xt3-0.1)))-2.1;

if(xp1 < xtp1 && xp2 < xtp2 && xp3 < xtp3){
	      xz1=test->vz[i1];
	      xz2=test->vz[i2];
	      xz3=test->vz[i3];
	      xvkey1=1;
	      xvkey2=1;
	      if(xz1 > xz3+3.0 || xz1 < xz3-3.0){
		 xvkey1=0;
	      }
	      if(xz2 > xz3+3.0 || xz2 < xz3-3.0){
		 xvkey2=0;
	      }	      
	      if(xz1==1 || xz2==1){

P4target.SetPxPyPzE(0.0,0.0,0.0,xmpr);
                 P4pho.SetPxPyPzE(0.,0.,egam,egam);
		 p_x=test->p[i1]*test->cx[i1];
		 p_y=test->p[i1]*test->cy[i1];
		 p_z=test->p[i1]*test->cz[i1];
		 e=sqrt(xmpi*xmpi+p_x*p_x+p_y*p_y+p_z*p_z);
                 P4pip1.SetPxPyPzE(p_x,p_y,p_z,e);
		 p_x=test->p[i2]*test->cx[i2];
		 p_y=test->p[i2]*test->cy[i2];
		 p_z=test->p[i2]*test->cz[i2];
		 e=sqrt(xmpi*xmpi+p_x*p_x+p_y*p_y+p_z*p_z);
                 P4pip2.SetPxPyPzE(p_x,p_y,p_z,e);
		 p_x=test->p[i3]*test->cx[i3];
		 p_y=test->p[i3]*test->cy[i3];
		 p_z=test->p[i3]*test->cz[i3];
		 e=sqrt(xmpi*xmpi+p_x*p_x+p_y*p_y+p_z*p_z);
                 P4pim.SetPxPyPzE(p_x,p_y,p_z,e);
		  
		 P4k01=P4pip1+P4pim;
		 P4k01=P4pip2+P4pim;
		 P4ksp=P4pip1+P4pip2+P4pim;
                 P4tot = P4target + P4pho;

xmk01=P4k01.M();
		 xmk02=P4k02.M();
		 xksp=P4ksp.M();
		 xmm=(P4tot-P4ksp).M();
		 ecm=sqrt(xmpr*xmpr+2*egam*xmpr);
		 beta=egam/(egam+xmpr);
		 gam=(egam+xmpr)/ecm;

pkpzcm=-gam*beta*P4ksp[3]+gam*P4ksp[2];
		 pkpcm=sqrt(P4ksp[0]*P4ksp[0]+P4ksp[1]*P4ksp[1]+pkpzcm*pkpzcm);
		 xang=pkpzcm/pkpcm;

mm=int(10*egam-15);
		 nn=int(10*(xang+1)/2.0*0.9);

XMK01->Fill(xmk01);
		 XMK02->Fill(xmk02);
		 if(fabs(xksp-xmksp)< cutksp){
		    if(fabs(xmk01-xmk0+2*cutk0)<cutk0||fabs(xmk01-xmk0+2*cutk0)<cutk0){
		       XKSPA->Fill(xksp); 
		       XMMA->Fill(xmm);
		       Lambda_A[nn][mm]->Fill(xmm);
		    }
		    if(fabs(xmk01-xmk0)<cutk0 ||fabs(xmk01-xmk0)<cutk0){
		       XKSPB->Fill(xksp);  
                       XMMB->Fill(xmm);
//		       Lambda_B[nn][mm]->Fill(xmm);
		    }
		    if(fabs(xmk01-xmk0-2*cutk0)<cutk0||fabs(xmk01-xmk0-2*cutk0)<cutk0){
		       XKSPC->Fill(xksp); 
		       XMMC->Fill(xmm);
//		       Lambda_C[nn][mm]->Fill(xmm);
		   }
		 }
	      }//vertex
           }//fiducial phi
	}//fiducial theta
      }//iselect
   }//event

TCanvas *c=new TCanvas("Histos","Histos",600,800);
   c->Divide(3,3);
   c->cd(1); TOF1->Draw();
   c->cd(2); TOF2->Draw();
   c->cd(3); TOF3->Draw();
   c->cd(4); TOF4->Draw();
   c->cd(5); TOF5->Draw();
   c->cd(6); TOF6->Draw();  
   c->cd(7); Lambda_A[5][23]->Draw();

Hlist->Write();
   outFile.Write(); // write to the output file
   outFile.Close(); // close the output file

}//end of main

//      TLorentzVector P4pim,P4pip1,P4pip2;
//      TLorentzVector P4target,P4pho,P4k0;
//      TLorentzVector P4tot,P4ksp;
//      Float_t p_x,p_y,p_z,e,pho4;

//      P4target.SetPxPyPzE(0.0,0.0,0.0,0.93828);
//      P4pho.SetPxPyPzE(0.,0.,pho4,pho4);
//      P4pim.SetPxPyPzE(p_x,p_y,p_z,e);
//      P4k0 = P4pip1+P4pim;
//      Float_t IM = P4k0.M();
//      P4ksp = P4k0+P4pip2;
//      P4tot = P4target + P4pho;
//      Float_t MMoffKs = (P4tot-P4ksp).M();

///////////////////////////////////////////////////////////////////////////////

float pg_corr(float p_in)
{
      Float_t x;Float_t pg_c;
      Float_t pc[3]={0.87854,0.42894E-01,0.85344};
      x=p_in/0.272;
      if(x>1) x=1;
      pg_c=pc[1]+pc[2]/pow(x,pc[3]);	
      return(pg_c);

}

int ibadpd(int isct, int ispd){

Int_t isect[11]={1,2,2,2,3,3,4,4,4,5,6};
   Int_t ipd[11]={33,8,17,27,18,19,8,13,14,23,30};
   Int_t ibad;

ibad=0;
   
   for(Int_t i=0; i<11; i++){
      if(isct == isect[i]){
         if(ispd == ipd[i]){
	    ibad=1;
         }
      }
   }
   return(ibad);
}

float xmm3(int i1,int i2,int i3,float xm1,float xm2,float xm3,float egam_in)
{

Float_t index[3]={-1};
      Float_t xm[3]={0};

index[0]=i1;
      index[1]=i2;
      index[2]=i3;
      xm[0] = xm1;
      xm[1] = xm2;
      xm[2] = xm3;
	
      Float_t xmm = -1.0;
      Float_t p0x = 0;
      Float_t p0y = 0;
      Float_t p0z = 0;
      Float_t ex  = 0;      
      Float_t e0  = 0;

for(Int_t i=0;i<3;i++){

Int_t ii=index[i];
     if(ii != -1){

p0x=p0x+(test->p[ii])*(test->cx[ii]);
       p0y=p0y+(test->p[ii])*(test->cy[ii]);
       p0z=p0z+(test->p[ii])*(test->cz[ii]);
       e0=sqrt((test->p[ii])*(test->p[ii])+xm[i]*xm[i]);		
       ex=ex+e0;

}
}     
     p0z=p0z-egam_in;

//  full momentum calculation

Float_t psq = p0x*p0x+p0y*p0y+p0z*p0z;			
      Float_t emiss = (egam_in + xmpr) - ex;
      Float_t x_x = emiss*emiss - psq;

xmm = sqrt( fabs(x_x) ); 
      
      if( x_x < 0.0 ) xmm = -xmm;

return(xmm);
}
float invm(int i1,int i2,int i3,float xm1,float xm2,float xm3)
{

Float_t index[3]={-1};
	Float_t xm[3]={0};

index[0]=i1;
	index[1]=i2;
	index[2]=i3;
	xm[0] = xm1;
	xm[1] = xm2;
	xm[2] = xm3;
	
      Float_t xmass = -1.0;
      Float_t p0x = 0;
      Float_t p0y = 0;
      Float_t p0z = 0;
      Float_t ex  = 0;      
      Float_t e0  = 0;

for(Int_t i=0;i<3;i++){

Int_t ii=index[i];
     if(ii != -1){

}
}

Float_t psq=p0x*p0x+p0y*p0y+p0z*p0z;			
      Float_t x_x = ex*ex - psq;

Float_t xmm = sqrt( fabs(x_x) ); 
      
      if( x_x > 0.0 ) xmass = sqrt(x_x);

return(xmass);
}

Example with Kinematic Fitting

int main(int __argc,char *__argv[]){
   
 char *outFileName = "k1k3fittest.root";
 extern int optind;   
 TROOT troot();  
 TFile outFile(outFileName,"recreate");  
 
  Float_t P_x[4],P_y[4],P_z[4],P_e[4],Vx[4],Vy[4],Vz[4],Cx[4],Cy[4],Cz[4];
 Int_t Q[4],Id[4];
 Float_t e_gamma,egamma;
  Float_t C11tber[4],C12tber[4],C13tber[4],C14tber[4],C15tber[4],C22tber[4],C23tber[4];
  Float_t C24tber[4],C25tber[4],C33tber[4],C34tber[4],C35tber[4],C44tber[4],C45tber[4],C55tber[4];

TLorentzVector P4pim1,P4pim2,P4pip1,P4pip2,P4kp,P4pr,P4target,P4pho,P4inc,P4allpi,P4mpr; 
 TLorentzVector Ppim1,Ppim2,Ppip1,Ppip2,Pkp,Ppho,Pinc,Pallpi; 
 TLorentzVector P4pho_K1,P4pim1_K1,P4pim2_K1,P4pip1_K1,P4pip2_K1,P4inc_K1,P4allpi_K1;
 TLorentzVector P4pho_K3,P4pim1_K3,P4pim2_K3,P4pip1_K3,P4pip2_K3,P4inc_K3,P4allpi_K3;
 TLorentzVector P4mpr_K3,P4mpr_K1;
 TLorentzVector P4k1,P4k2,P4k3,P4k4;
 TLorentzVector Pk1,Pk2,Pk3,Pk4;
 TVector3 P3pip1, P3pip2, P3pim1, P3pim2, C3pip1,C3pip2,C3pim1,C3pim2;
 TVector3 V3pip1,V3pip2,V3pim1,V3pim2,K1track,K2track,K3track,K4track;
 TVector3 F1track,F2track,P3k1,P3k2,P3k3,P3k4;
 Float_t xmm, xmm1,xmm2,xk1m,xk2m,xk3m,xk4m,xm4pi_k12, xm4pi_k34;    
 Float_t xmmk1,xmmk2,xmmk3,xmmk4;
 Float_t PLk1,PLk2,PLk3,PLk4,Betak1,Betak2,Betak3,Betak4,Ltimek1,Ltimek2,Ltimek3,Ltimek4;
 Double_t xm2k1p,xm2k2p,xm2k3p,xm2k4p,tvar_k12,tvar_k34;
 Float_t distp1m1,distp2m2,distp1m2,distp2m1,distk1k2,distk3k4;
 Float_t egam, costh_k12,costh_k34,egam_K1,egam_K3;
 Float_t xright,xleft;
// DOCA K1,K2,K3,K4,F1,F2;  //struct DOCA { Float_t distance; TVector3 vertex;} defined in doca.c   
//DOCA doca(TVector3,TVector3,TVector3,TVector3);
 float coscm(TLorentzVector,TLorentzVector,TLorentzVector,TLorentzVector, Float_t);

const Float_t kkincut =0.02;
 const Float_t k0cut = 0.008386;   //0.009614;
 const Float_t sd = 2.0;
 const Float_t xmks=0.497614;
 const Float_t xmpr=0.938272;

//  Various Histograms
 TH1F *prob = new TH1F("prob","prob",100,-0.01,1.1);
 TH1F *chi2 = new TH1F("chi2","#chi^{2}",100,0.0,5.0);
 TH1F *prob_K1 = new TH1F("prob_K1","prob",100,-0.01,1.1);
 TH1F *chi2_K1 = new TH1F("chi2_K1","#chi^{2}",100,0.0,5.0);
 TH1F *prob_K3 = new TH1F("prob_K3","prob",100,-0.01,1.1);
 TH1F *chi2_K3 = new TH1F("chi2_K3","#chi^{2}",100,0.0,5.0);
 TH1F *new_mm = new TH1F("fit_mm","fit_mm",100,0.6,1.4);
 
 TH1F *Pull_Egam = new TH1F("Pull_Egam","#gamma E-pull",100,-6.0,6.0);
 TH1F *Pull_pip1p = new TH1F("Pull_pip1p","#pi^{+}_{fast} p-Pull",100,-6.0,6.0);
 TH1F *Pull_pip1lam = new TH1F("Pull_pip1lam","#pi^{+}_{fast} #lambda-pull",100,-6.0,6.0);
 TH1F *Pull_pip1phi = new TH1F("Pull_pip1phi","#pi^{+}_{fast} #phi-pull",100,-6.0,6.0);
 TH1F *Pull_pip2p = new TH1F("Pull_pip2p","#pi^{+}_{slow} p-pull",100,-6.0,6.0);
 TH1F *Pull_pip2lam = new TH1F("Pull_pip2lam","#pi^{+}_{slow} #lambda-pull",100,-6.0,6.0);
 TH1F *Pull_pip2phi= new TH1F("Pull_pip2phi","#pi^{+}_{slow} #phi-pull",100,-6.0,6.0);
 TH1F *Pull_pim1p = new TH1F("Pull_pim1p","#pi^{-}_{fast} p-pull",100,-6.0,6.0);
 TH1F *Pull_pim1lam = new TH1F("Pull_pim1lam","#pi^{-}_{fast} #lambda-pull",100,-6.0,6.0);
 TH1F *Pull_pim1phi = new TH1F("Pull_pim1phi","#pi^{-}_{fast} #phi-pull",100,-6.0,6.0);
 TH1F *Pull_pim2p = new TH1F("Pull_pim2p","#pi^{-}_{slow} p-pull",100,-6.0,6.0);
 TH1F *Pull_pim2lam = new TH1F("Pull_pim2lam","#pi^{-}_{slow} #lambda-pull",100,-6.0,6.0);
 TH1F *Pull_pim2phi = new TH1F("Pull_pim2phi","#pi^{-}_{slow} #phi-pull",100,-6.0,6.0);
 
 for(int n_arg = optind; n_arg < __argc; n_arg++){   
   TFile inFile(__argv[n_arg]); // open the input file         
   
   if( TTree *p0 = (TTree*)inFile.Get("p0")){

//p0->SetBranchAddress("nrun",&nrun); 
    p0->SetBranchAddress("Id", Id);
    p0->SetBranchAddress("e_gamma", &e_gamma);
//  p0->SetBranchAddress("Vertex_time", &Vertex_time);
//  p0->SetBranchAddress("xmm",&xmm);
    p0->SetBranchAddress("P_x",P_x);
    p0->SetBranchAddress("P_y",P_y);
    p0->SetBranchAddress("P_z",P_z);
    p0->SetBranchAddress("P_e",P_e);
    p0->SetBranchAddress("Vx",Vx);
    p0->SetBranchAddress("Vy",Vy);
    p0->SetBranchAddress("Vz",Vz);
    p0->SetBranchAddress("Cx",Cx);
    p0->SetBranchAddress("Cy",Cy);
    p0->SetBranchAddress("Cz",Cz);
    p0->SetBranchAddress("C11tber",C11tber);
    p0->SetBranchAddress("C12tber",C12tber);
    p0->SetBranchAddress("C13tber",C13tber);
    p0->SetBranchAddress("C14tber",C14tber);
    p0->SetBranchAddress("C15tber",C15tber);
    p0->SetBranchAddress("C22tber",C22tber);
    p0->SetBranchAddress("C23tber",C23tber);
    p0->SetBranchAddress("C24tber",C24tber);
    p0->SetBranchAddress("C25tber",C25tber);
    p0->SetBranchAddress("C33tber",C33tber);
    p0->SetBranchAddress("C34tber",C34tber);
    p0->SetBranchAddress("C35tber",C35tber);
    p0->SetBranchAddress("C44tber",C44tber);
    p0->SetBranchAddress("C45tber",C45tber);
    p0->SetBranchAddress("C55tber",C55tber);

Int_t nentries = (Int_t)p0->GetEntries();  
     for (Int_t j=0;j<nentries;j++) { p0->GetEntry(j);
        P4target.SetPxPyPzE(0.0,0.0,0.0,xmpr);
        egamma = e_gamma;
        Ppho.SetPxPyPzE(0.,0.,egamma,egamma);
     	Pinc =  P4target + Ppho;

Ppip1.SetPxPyPzE(P_x[0],P_y[0],P_z[0],P_e[0]); 
        Ppip2.SetPxPyPzE(P_x[1],P_y[1],P_z[1],P_e[1]);
        Ppim1.SetPxPyPzE(P_x[2],P_y[2],P_z[2],P_e[2]); 
        Ppim2.SetPxPyPzE(P_x[3],P_y[3],P_z[3],P_e[3]);

Pallpi = Ppip1+Ppip2+Ppim1+Ppim2;

P4mpr = Pinc-Pallpi;
        xmm = fabs(P4mpr.M());
  
        V3pip1.SetXYZ(Vx[0],Vy[0],Vz[0]);       V3pip2.SetXYZ(Vx[1],Vy[1],Vz[1]);
        V3pim1.SetXYZ(Vx[2],Vy[2],Vz[2]);      V3pim2.SetXYZ(Vx[3],Vy[3],Vz[3]);

C3pip1.SetXYZ(Cx[0],Cy[0],Cz[0]);       C3pip2.SetXYZ(Cx[1],Cy[1],Cz[1]);
        C3pim1.SetXYZ(Cx[2],Cy[2],Cz[2]);      C3pim2.SetXYZ(Cx[3],Cy[3],Cz[3]);

Pk1 = Ppip1+Ppim1; 
        Pk2 = Ppip2 + Ppim2; 
        Pk3 = Ppip1 + Ppim2; 
        Pk4 = Ppip2 + Ppim1;

xk1m = Pk1.M(); //Invariant mass of K1
        xk2m = Pk2.M(); //Invariant mass of K2
        xk3m = Pk3.M(); //Invariant mass of K3
        xk4m = Pk4.M(); //Invariant mass of K4

//--------------------------------------------------------Kinematic Fitting---------------------------------------------------
//          if(fabs(xmm-xmpr)<0.05) {	 
//   	  if(  ( (fabs(xk1m-xmks)<kkincut) && (fabs(xk2m-xmks)<kkincut) )  || ( (fabs(xk3m-xmks)<kkincut) && (fabs(xk4m-xmks)<kkincut) )  )  {

const int nparts = 4;  //number of particles in fit
	 	  
	  vector<TLorentzVector> Vin4(nparts);  //set of 4-Vectors use in the fit
	  std::vector<TVector3> Vert(nparts);  //set of verties 
	  std::vector<string> p_names(nparts);  // string names
	 
	  /*This next bool sets the multiscattering and energy spread correction
	  algorithm on (true) or off (false).  If off is used the CMU averaging
	  technique is employed that reads files for particular particle types.
	  If tuned correctly the averaging method can work just as good or maybe better
	  for a reaction that has all final decay products starting at the primary
	  interaction vertex.  Anything that has a secondary decay vertex can lead
	  to inconsistent covariance and flag (true) is required.*/
	  
	  
	  bool Multi = false; 
	  	  
	  /*This next bool sets the multiscattering and energy spread correction
	  algorithm on (true) or off (false).  If off is used the CMU averaging
	  technique is employed that reads files form the ./parms directory for
	  particular particle types.  If tuned correctly the averaging method can 
	  work just as good or maybe better for a reaction that has all final decay 
	  products starting at the primary interaction vertex.  Anything that decays 
	  from a secondary decay vertex can lead to inconsistent covariance and flag 
	  (true) is required.  Option false requires that you have tuned parameters
	  for you particle species.  These parms are called in KinLine.C and needs
	  to point to your copy of parms: */

//  String fileName_base = "/home/cshloka/Makeshift/package/parms/";	  
	  bool MC = false;
	   
	  //This is for extra constraints not used in this example
	  
	  std::vector<bool> set(nparts);
	  set[0] = false;
	  set[1] = false;			      
 	  set[2] = false;
	  set[3] = false;

Double_t m_targ = 0.938272;    //target mass
	  const int dim = 3*nparts+1;    //dimension of covariance matrix
	  TMatrixD covMatrix(dim,dim); //Matrix you will use
	  TMatrixD covTrack(dim,dim);  // Matrix from CLAS TBER 
          Float_t p_mag[4];
	  
	  // detected 4-Vectors
	  Vin4[0] = Ppip1;   p_mag[0] = Ppip1.P();
	  Vin4[1] = Ppip2;   p_mag[1] = Ppip2.P();
	  Vin4[2] = Ppim1;  p_mag[2] = Ppim1.P();
	  Vin4[3] = Ppim1;  p_mag[3] = Ppim2.P();

// vertex information or best guess
	  // this is used in Multiple Scattering and energy spread
	  Vert[0] = V3pip1;
	  Vert[1] = V3pip2;
	  Vert[2] = V3pim1;
	  Vert[3] = V3pim2;

// for particle names use pi-, pi+,e+,e-,k+,k-,(for gamma see HowToPhotons.txt)
	  //These must be filled in the array in the same order
	  //as the covariance matrix for each particle	  
	  p_names[0] = "pi+";
	  p_names[1] = "pi+";
  	  p_names[2] = "pi-";
	  p_names[3] = "pi-";
 
 /* To fill covTrack you need the first three
	  components of the TBER information c11,c22,c33,c12,c13,c23,c32

covTrack(0,0) don't fill */
   	  for (Int_t i =0; i<4; i++) {
   	  covTrack(3*i+1,3*i+1) = C11tber[i]*pow(p_mag[i],4); 
   	  covTrack(3*i+2,3*i+2) = C22tber[i];
   	  covTrack(3*i+3,3*i+3) = C33tber[i];
   	  covTrack(3*i+1,3*i+2) =0;// -C12tber[i]*Q[i]*p_mag[i]*p_mag[i];
   	  covTrack(3*i+2,3*i+1) =0;// -C12tber[i]*Q[i]*p_mag[i]*p_mag[i];
   	  covTrack(3*i+1,3*i+3) =0;// -C13tber[i]*Q[i]*p_mag[i]*p_mag[i];
   	  covTrack(3*i+3,3*i+1) =0;// -C13tber[i]*Q[i]*p_mag[i]*p_mag[i];
   	  covTrack(3*i+2,3*i+3) = C23tber[i];
   	  covTrack(3*i+3,3*i+2) = C23tber[i]; 
cout<<Q[i]<<endl;	  
	  }

// Get The full covariance for your fit
	  string experiment = "g11a";       //options  g10a,g11a,g12,g13a
	  covMatrix= CorrectCLAS_V(covTrack,p_names,Vin4,Vert,Multi,MC,experiment);
//covTrack.Print();
	  // Declare your fit	  
	  Kstream FirstFit;
	  FirstFit.StringNames(p_names);
	  FirstFit.FitInput(egamma,Vin4,covMatrix,m_targ);	  
	  // Declare the fit Hypothesis (use particle names a listed)
	  FirstFit.Fit(); 	  // missing pi+

prob->Fill(FirstFit.Prob());  //x:0-1
	 chi2->Fill(FirstFit.Chi2());  //x:0-20
  if(FirstFit.Prob() > 0.01){
	  //Fill histos with stretch info

Pull_Egam->Fill(FirstFit.GetPull(0));
	  Pull_pip1p->Fill(FirstFit.GetPull(1));	  
	  Pull_pip1lam->Fill(FirstFit.GetPull(2));	  
	  Pull_pip1phi->Fill(FirstFit.GetPull(3));
	  Pull_pip2p->Fill(FirstFit.GetPull(4));	  
	  Pull_pip2lam->Fill(FirstFit.GetPull(5));	  
	  Pull_pip2phi->Fill(FirstFit.GetPull(6));
	  Pull_pim1p->Fill(FirstFit.GetPull(7));	  
	  Pull_pim1lam->Fill(FirstFit.GetPull(8));	  
	  Pull_pim1phi->Fill(FirstFit.GetPull(9));
	  Pull_pim2p->Fill(FirstFit.GetPull(10));	  
	  Pull_pim2lam->Fill(FirstFit.GetPull(11));	  
	  Pull_pim2phi->Fill(FirstFit.GetPull(12));

}

//      } //loose cuts on K masses
//      }  //cut on missing mass for KF

}//for entries
   }//if TTree *p0 = (TTree*)

else {
       cout << "File has no TTree named p0 or does not exist!!" << endl;
   }  
 }//n_arg

outFile.Write(); // write to the output file
 outFile.Close(); // close the output file
        
//end of main
}