[u/mrichter/AliRoot.git] / MUON / MUONTriggerEfficiencyPt.C

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */


// Macro to produce trigger single muon efficiency versus pt plots for the 
// 2 pt cuts. 
// see full description in the README file
// Author: Fabien Guerin (LPC)

// ROOT includes
#include "TBranch.h"
#include "TClonesArray.h"
#include "TFile.h"
#include "TF1.h"
#include "TH1.h"
#include "TH2.h"
#include "TParticle.h"
#include "TTree.h"
#include "TMath.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TLegend.h"
#include "Riostream.h"

// STEER includes
#include "AliRun.h"
#include "AliRunLoader.h"
#include "AliHeader.h"
#include "AliLoader.h"
#include "AliStack.h"

// MUON includes
#include "AliMUONSimData.h"
#include "AliMUONRecData.h"
#include "AliMUONHit.h"
#include "AliMUONDigit.h"
#include "AliMUONGlobalTrigger.h"
#include "AliMUONLocalTrigger.h"
#include "AliMUONTrack.h"

Double_t fitArc(Double_t *x,Double_t *par) 
{
  Double_t h1=par[1]*(x[0]+par[2]);
  return par[0]*TMath::ATan(TMath::Power(h1,6))+par[3];
}

Double_t fitArch(Double_t *x,Double_t *par) 
{
  Double_t h0=2*par[0]/(TMath::Pi());
  Double_t h1=par[1]*x[0]+par[2]*x[0]*x[0];
  return h0*TMath::TanH(h1)+par[3];
}

void MUONTriggerEfficiencyPt(char filenameSim[10]="galice.root",  
                             char filenameRec[10]="galice.root",
                             Bool_t readFromRP = 0)
{

// define style
    TStyle *st1 = new TStyle("st1","My STYLE");
    
    st1->SetOptStat(0);
    st1->SetOptFit(111);  
    st1->SetFrameFillColor(10);
    st1->SetCanvasColor(10);
    st1->SetFillColor(10);
    st1->SetTitleBorderSize(0);
    st1->SetTitleOffset(1.2,"XY");
    st1->SetTitleSize(0.05,"XY"); 
    st1->SetLabelSize(0.045,"XY");  
    st1->SetLabelFont(22,"XY");  
    st1->SetTitleFont(22,"XY");     
    st1->SetOptTitle(0);   
    st1->SetStatColor(10);
    st1->SetStatFont(62);     
    st1->SetFitFormat("4.4g");
    st1->SetPadLeftMargin(0.15);
    st1->SetPadRightMargin(0.06); 
    st1->SetPadTopMargin(0.06);
    st1->SetPadBottomMargin(0.15); 
    st1->cd();
    
//    gROOT->ForceStyle();
    //TGaxis::SetMaxDigits(3);
    
// beginning of macro    
    TH1F *ptlpt    = new TH1F("ptlpt","",50,0.15,5.); 
    TH1F *pthpt    = new TH1F("pthpt","",50,0.15,5.);
    TH1F *ptcoinc  = new TH1F("ptcoinc","",50,0.15,5.);  
        
    TParticle *particle;
    AliStack* stack; 
    
    Int_t nevents;
    Double_t coincmuon=0;
    Double_t lptmuon=0;
    Double_t hptmuon=0;
    Double_t ptmu=0;

// output file
    char digitdat[100];
    char currentfile[100];
    sprintf(digitdat,"MUONTriggerEfficiencyPt.out");  
    FILE *fdat=fopen(digitdat,"w");
    
// Initialise AliRoot
// Creating Run Loader and openning file containing Hits
    AliRunLoader * RunLoaderSim = AliRunLoader::Open(filenameSim,"MUONFolderSim","READ");
    
    if (RunLoaderSim ==0x0) {
	printf(">>> Error : Error Opening %s file \n",currentfile);
	return;
    }
    
    AliRunLoader * RunLoaderRec = AliRunLoader::Open(filenameRec,"MUONFolder","READ");
    
    if (RunLoaderRec ==0x0) {
	printf(">>> Error : Error Opening %s file \n",currentfile);
	return;
    }
    
    nevents = RunLoaderSim->GetNumberOfEvents();
    
    AliLoader * MUONLoaderSim = RunLoaderSim->GetLoader("MUONLoader");
    AliLoader * MUONLoaderRec = RunLoaderRec->GetLoader("MUONLoader");
    
    if (!readFromRP) {
	cout << " reading from digits \n";
	MUONLoaderSim->LoadDigits("READ");
    } else {
	cout << " reading from RecPoints \n";
	MUONLoaderRec->LoadRecPoints("READ");
   }
    MUONLoaderSim->LoadHits("READ");
    RunLoaderSim->LoadKinematics("READ");


// Creating MUON data containers
    AliMUONSimData muondataSim(MUONLoaderSim,"MUON","MUON");
    AliMUONRecData muondataRec(MUONLoaderRec,"MUON","MUON");

    TClonesArray * globalTrigger;
    AliMUONGlobalTrigger * gloTrg;

    for (Int_t ievent=0; ievent<nevents; ievent++) {  // Event loop
	
	RunLoaderSim->GetEvent(ievent);
	RunLoaderRec->GetEvent(ievent);
	
	if (ievent%500==0) printf("ievent = %d \n",ievent);

// kine
	Int_t iparticle, nparticles;
	stack = RunLoaderSim->Stack();
	nparticles = (Int_t) stack->GetNtrack();        
	for (iparticle=0; iparticle<nparticles; iparticle++) {             
            particle = stack->Particle(iparticle);   
            Float_t pt = particle->Pt();	    
            Int_t pdgcode = particle->GetPdgCode();                        
	    if (pdgcode==-13 || pdgcode==13) ptmu = pt;             
        }

// trigger 
	if (!readFromRP) {
	    muondataSim.SetTreeAddress("D,GLT"); 
	    muondataSim.GetTriggerD();
	    globalTrigger = muondataSim.GlobalTrigger();
	} else {    
	    muondataRec.SetTreeAddress("RC,TC"); 
	    muondataRec.GetTrigger();
	    globalTrigger = muondataRec.GlobalTrigger();

	}

        Int_t nglobals = (Int_t) globalTrigger->GetEntriesFast(); // should be 1 
        for (Int_t iglobal=0; iglobal<nglobals; iglobal++) { // Global Trigger
	    gloTrg = static_cast<AliMUONGlobalTrigger*>(globalTrigger->At(iglobal));
	    
	    if (gloTrg->SingleLpt()>=1 ) {
		lptmuon++;
		ptlpt->Fill(ptmu);
	    }       
	    if (gloTrg->SingleHpt()>=1 ) {
		hptmuon++;
		pthpt->Fill(ptmu);
	    }
        } // end of loop on Global Trigger      
        muondataSim.ResetTrigger();  
        muondataRec.ResetTrigger();  

// Hits
	muondataSim.SetTreeAddress("H");    

        Int_t itrack, ntracks, NbHits[4];
        Int_t SumNbHits;

        for (Int_t j=0; j<4; j++) NbHits[j]=0;
 
	ntracks = (Int_t) muondataSim.GetNtracks();
   
        for (itrack=0; itrack<ntracks; itrack++) { // track loop
	    muondataSim.GetTrack(itrack);

         Int_t ihit, nhits;
	 nhits = (Int_t) muondataSim.Hits()->GetEntriesFast();
         AliMUONHit* mHit;
    
          for (ihit=0; ihit<nhits; ihit++) {
	      mHit = static_cast<AliMUONHit*>(muondataSim.Hits()->At(ihit));
            Int_t Nch        = mHit->Chamber(); 
            Int_t hittrack   = mHit->Track();
            Float_t IdPart   = mHit->Particle();	    

            for (Int_t j=0;j<4;j++) {
              Int_t kch=11+j;
	      if (hittrack==0) {
                if (Nch==kch && (IdPart==-13 || IdPart==13) && NbHits[j]==0) NbHits[j]++; 
		}               
            }    
          }
	  muondataSim.ResetHits();
        } // end track loop

// 3/4 coincidence 
        SumNbHits=NbHits[0]+NbHits[1]+NbHits[2]+NbHits[3];

        if (SumNbHits==3 || SumNbHits==4) {
            coincmuon++;
            ptcoinc->Fill(ptmu);
        }            
                
    } // end loop on event
    
    if (coincmuon==0) {
	cout << " >>> <E> coincmuon = 0 after event loop " << "\n";
	cout << " >>> this probably means that input does not contain one (and only one) " << "\n";
        cout << " >>> muon track per event as it should " << "\n";
	cout << " >>> see README for further information " << "\n";
	cout << " >>> exiting now ! " << "\n";
	return;
    }

      
       MUONLoaderSim->UnloadHits();
       if (!readFromRP) {
	   MUONLoaderSim->UnloadDigits();  
       } else {    
	   MUONLoaderRec->UnloadRecPoints();
       }
       RunLoaderSim->UnloadKinematics();
     
      fprintf(fdat,"\n");    
      fprintf(fdat,"\n");
      fprintf(fdat," Number of events = %d \n",nevents);  
      fprintf(fdat," Number of events with 3/4 coinc = %d \n",(Int_t)coincmuon);
      fprintf(fdat," Nomber of dimuons with 3/4 coinc Lpt cut = %d \n",(Int_t)lptmuon); 
      fprintf(fdat," Number of dimuons with 3/4 coinc Hpt cut = %d \n",(Int_t)hptmuon);  
      fprintf(fdat,"\n"); 
      
      Double_t efficiency,error;

      efficiency=lptmuon/coincmuon;  
      error=efficiency*TMath::Sqrt((lptmuon+coincmuon)/(lptmuon*coincmuon));  
      fprintf(fdat," Efficiency Lpt cut = %4.4f +/- %4.4f\n",efficiency,error);
      
      efficiency=hptmuon/coincmuon; 
      error=efficiency*TMath::Sqrt((hptmuon+coincmuon)/(hptmuon*coincmuon)); 
      fprintf(fdat," Efficiency Hpt cut = %4.4f +/- %4.4f\n",efficiency,error);
      
      fclose(fdat);
      
      Double_t x1,x2,xval,xerr;
      
      for (Int_t i=0;i<50;i++) {
	  x1=ptcoinc->GetBinContent(i+1);
	  
	  x2=ptlpt->GetBinContent(i+1);
	  if (x1!=0 && x2!=0) {
	      xval=x2/x1;
	      xerr=xval*TMath::Sqrt((x1+x2)/x1*x2);   
	      ptlpt->SetBinContent(i+1,xval);
	      ptlpt->SetBinError(i+1,0);
	  } else {
	      ptlpt->SetBinContent(i+1,0.);
	      ptlpt->SetBinError(i+1,0.);     
	  } 
	  
	  x2=pthpt->GetBinContent(i+1);
	  if (x1!=0 && x2!=0) {
	      xval=x2/x1;
	      xerr=xval*TMath::Sqrt((x1+x2)/x1*x2);   
	      pthpt->SetBinContent(i+1,xval);
	      pthpt->SetBinError(i+1,0);
	  } else {
	      pthpt->SetBinContent(i+1,0.);
	      pthpt->SetBinError(i+1,0.);     
	  }              
      }      
      
      TF1 *fitlpt = new TF1("fitlpt",fitArc,0.,5.,4); 
      TF1 *fithpt = new TF1("fithpt",fitArc,0.,5.,4);      
      
      TCanvas *c1 = new TCanvas("c1","Trigger efficiency vs pt and y for single muon",200,0,900,400); 
      c1->Divide(2,1);      
      
      c1->cd(1);
      ptlpt->SetMinimum(0.);
      ptlpt->SetMaximum(1.05);   
      ptlpt->SetTitle("");
      ptlpt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
      ptlpt->GetYaxis()->SetTitle("Efficiency");     
      //ptlpt->GetXaxis()->SetRange(3);
      ptlpt->Draw(""); 
      fitlpt->SetParameters(0.602,0.774,0.273,0.048);  
      fitlpt->SetLineColor(2);
      fitlpt->SetLineWidth(2);
      fitlpt->Draw("SAME");    
      TLegend * leg = new TLegend(0.5,0.38,0.85,0.53); 
      leg = new TLegend(0.5,0.38,0.85,0.53); 
      leg->SetBorderSize(0);
      leg->SetFillColor(0);
      leg->SetTextSize(0.05);
      leg->SetTextFont(22);
      leg->SetHeader("Lpt trigger pt cut");
      leg->AddEntry(fitlpt,"Ref.","l");
      leg->AddEntry(ptlpt,"New","l");                 
      leg->Draw("SAME");
      
      c1->cd(2);
      pthpt->SetMinimum(0.);
      pthpt->SetMaximum(1.05);   
      pthpt->SetTitle("");
      pthpt->GetXaxis()->SetTitle("P_{T} (GeV/c)");  
      pthpt->GetYaxis()->SetTitle("Efficiency");     
      //pthpt->GetXaxis()->SetRange(3);
      pthpt->Draw(""); 
      fithpt->SetParameters(0.627,0.393,0.855,0.0081); 
      fithpt->SetLineColor(2);
      fithpt->SetLineWidth(2);
      fithpt->Draw("SAME");    
      leg = new TLegend(0.5,0.38,0.85,0.53); 
      leg->SetBorderSize(0);
      leg->SetFillColor(0);
      leg->SetTextSize(0.05);
      leg->SetTextFont(22);
      leg->SetHeader("Hpt trigger pt cut");
      leg->AddEntry(fithpt,"Ref.","l");
      leg->AddEntry(pthpt,"New","l");                 
      leg->Draw("SAME");
      
      c1->SaveAs("MUONTriggerEfficiencyPt.gif");
      c1->SaveAs("MUONTriggerEfficiencyPt.eps");      
}
Commit	Line	Data
5707a463	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	/* $Id$ */
	17
	18
07d83718	19	// Macro to produce trigger single muon efficiency versus pt plots for the
	20	// 2 pt cuts.
	21	// see full description in the README file
	22	// Author: Fabien Guerin (LPC)
5707a463	23
	24	// ROOT includes
	25	#include "TBranch.h"
	26	#include "TClonesArray.h"
	27	#include "TFile.h"
	28	#include "TF1.h"
	29	#include "TH1.h"
	30	#include "TH2.h"
	31	#include "TParticle.h"
	32	#include "TTree.h"
	33	#include "TMath.h"
	34	#include "TStyle.h"
	35	#include "TCanvas.h"
	36	#include "TLegend.h"
b19dcb40	37	#include "Riostream.h"
5707a463	38
	39	// STEER includes
	40	#include "AliRun.h"
5707a463	41	#include "AliRunLoader.h"
	42	#include "AliHeader.h"
	43	#include "AliLoader.h"
	44	#include "AliStack.h"
5707a463	45
5707a463	46	// MUON includes
b19dcb40	47	#include "AliMUONSimData.h"
b19dcb40	48	#include "AliMUONRecData.h"
5707a463	49	#include "AliMUONHit.h"
5707a463	50	#include "AliMUONDigit.h"
5707a463	51	#include "AliMUONGlobalTrigger.h"
	52	#include "AliMUONLocalTrigger.h"
	53	#include "AliMUONTrack.h"
	54
	55	Double_t fitArc(Double_t x,Double_t par)
	56	{
	57	Double_t h1=par[1]*(x[0]+par[2]);
	58	return par[0]*TMath::ATan(TMath::Power(h1,6))+par[3];
	59	}
	60
	61	Double_t fitArch(Double_t x,Double_t par)
	62	{
	63	Double_t h0=2*par[0]/(TMath::Pi());
	64	Double_t h1=par[1]x[0]+par[2]x[0]*x[0];
	65	return h0*TMath::TanH(h1)+par[3];
	66	}
	67
b19dcb40	68	void MUONTriggerEfficiencyPt(char filenameSim[10]="galice.root",
	69	char filenameRec[10]="galice.root",
	70	Bool_t readFromRP = 0)
5707a463	71	{
	72
	73	// define style
	74	TStyle *st1 = new TStyle("st1","My STYLE");
	75
	76	st1->SetOptStat(0);
	77	st1->SetOptFit(111);
	78	st1->SetFrameFillColor(10);
	79	st1->SetCanvasColor(10);
	80	st1->SetFillColor(10);
	81	st1->SetTitleBorderSize(0);
	82	st1->SetTitleOffset(1.2,"XY");
	83	st1->SetTitleSize(0.05,"XY");
	84	st1->SetLabelSize(0.045,"XY");
	85	st1->SetLabelFont(22,"XY");
	86	st1->SetTitleFont(22,"XY");
	87	st1->SetOptTitle(0);
	88	st1->SetStatColor(10);
	89	st1->SetStatFont(62);
	90	st1->SetFitFormat("4.4g");
	91	st1->SetPadLeftMargin(0.15);
	92	st1->SetPadRightMargin(0.06);
	93	st1->SetPadTopMargin(0.06);
	94	st1->SetPadBottomMargin(0.15);
	95	st1->cd();
	96
07d83718	97	// gROOT->ForceStyle();
5707a463	98	//TGaxis::SetMaxDigits(3);
	99
	100	// beginning of macro
5707a463	101	TH1F *ptlpt = new TH1F("ptlpt","",50,0.15,5.);
	102	TH1F *pthpt = new TH1F("pthpt","",50,0.15,5.);
	103	TH1F *ptcoinc = new TH1F("ptcoinc","",50,0.15,5.);
	104
	105	TParticle *particle;
	106	AliStack* stack;
	107
b19dcb40	108	Int_t nevents;
5707a463	109	Double_t coincmuon=0;
5707a463	110	Double_t lptmuon=0;
5707a463	111	Double_t hptmuon=0;
522b4724	112	Double_t ptmu=0;
5707a463	113
	114	// output file
	115	char digitdat[100];
b19dcb40	116	char currentfile[100];
5707a463	117	sprintf(digitdat,"MUONTriggerEfficiencyPt.out");
5707a463	118	FILE *fdat=fopen(digitdat,"w");
b19dcb40	119
b19dcb40	120	// Initialise AliRoot
5707a463	121	// Creating Run Loader and openning file containing Hits
b19dcb40	122	AliRunLoader * RunLoaderSim = AliRunLoader::Open(filenameSim,"MUONFolderSim","READ");
	123
	124	if (RunLoaderSim ==0x0) {
	125	printf(">>> Error : Error Opening %s file \n",currentfile);
	126	return;
	127	}
5707a463	128
b19dcb40	129	AliRunLoader * RunLoaderRec = AliRunLoader::Open(filenameRec,"MUONFolder","READ");
	130
	131	if (RunLoaderRec ==0x0) {
	132	printf(">>> Error : Error Opening %s file \n",currentfile);
5707a463	133	return;
	134	}
	135
b19dcb40	136	nevents = RunLoaderSim->GetNumberOfEvents();
5707a463	137
b19dcb40	138	AliLoader * MUONLoaderSim = RunLoaderSim->GetLoader("MUONLoader");
b19dcb40	139	AliLoader * MUONLoaderRec = RunLoaderRec->GetLoader("MUONLoader");
5707a463	140
b19dcb40	141	if (!readFromRP) {
	142	cout << " reading from digits \n";
	143	MUONLoaderSim->LoadDigits("READ");
	144	} else {
	145	cout << " reading from RecPoints \n";
	146	MUONLoaderRec->LoadRecPoints("READ");
	147	}
	148	MUONLoaderSim->LoadHits("READ");
	149	RunLoaderSim->LoadKinematics("READ");
	150
	151
	152	// Creating MUON data containers
	153	AliMUONSimData muondataSim(MUONLoaderSim,"MUON","MUON");
	154	AliMUONRecData muondataRec(MUONLoaderRec,"MUON","MUON");
	155
5707a463	156	TClonesArray * globalTrigger;
	157	AliMUONGlobalTrigger * gloTrg;
	158
5707a463	159	for (Int_t ievent=0; ievent<nevents; ievent++) { // Event loop
5707a463	160
b19dcb40	161	RunLoaderSim->GetEvent(ievent);
b19dcb40	162	RunLoaderRec->GetEvent(ievent);
5707a463	163
	164	if (ievent%500==0) printf("ievent = %d \n",ievent);
	165
5707a463	166	// kine
5707a463	167	Int_t iparticle, nparticles;
b19dcb40	168	stack = RunLoaderSim->Stack();
5707a463	169	nparticles = (Int_t) stack->GetNtrack();
	170	for (iparticle=0; iparticle<nparticles; iparticle++) {
	171	particle = stack->Particle(iparticle);
	172	Float_t pt = particle->Pt();
	173	Int_t pdgcode = particle->GetPdgCode();
	174	if (pdgcode==-13 \|\| pdgcode==13) ptmu = pt;
	175	}
	176
	177	// trigger
07d83718	178	if (!readFromRP) {
b19dcb40	179	muondataSim.SetTreeAddress("D,GLT");
	180	muondataSim.GetTriggerD();
	181	globalTrigger = muondataSim.GlobalTrigger();
07d83718	182	} else {
b19dcb40	183	muondataRec.SetTreeAddress("RC,TC");
	184	muondataRec.GetTrigger();
	185	globalTrigger = muondataRec.GlobalTrigger();
	186
07d83718	187	}
5707a463	188
	189	Int_t nglobals = (Int_t) globalTrigger->GetEntriesFast(); // should be 1
	190	for (Int_t iglobal=0; iglobal<nglobals; iglobal++) { // Global Trigger
	191	gloTrg = static_cast<AliMUONGlobalTrigger*>(globalTrigger->At(iglobal));
	192
522b4724	193	if (gloTrg->SingleLpt()>=1 ) {
5707a463	194	lptmuon++;
	195	ptlpt->Fill(ptmu);
	196	}
522b4724	197	if (gloTrg->SingleHpt()>=1 ) {
5707a463	198	hptmuon++;
	199	pthpt->Fill(ptmu);
	200	}
	201	} // end of loop on Global Trigger
b19dcb40	202	muondataSim.ResetTrigger();
b19dcb40	203	muondataRec.ResetTrigger();
5707a463	204
5707a463	205	// Hits
b19dcb40	206	muondataSim.SetTreeAddress("H");
b19dcb40	207
5707a463	208	Int_t itrack, ntracks, NbHits[4];
5707a463	209	Int_t SumNbHits;
522b4724	210
5707a463	211	for (Int_t j=0; j<4; j++) NbHits[j]=0;
5707a463	212
b19dcb40	213	ntracks = (Int_t) muondataSim.GetNtracks();
b19dcb40	214
5707a463	215	for (itrack=0; itrack<ntracks; itrack++) { // track loop
b19dcb40	216	muondataSim.GetTrack(itrack);
5707a463	217
5707a463	218	Int_t ihit, nhits;
b19dcb40	219	nhits = (Int_t) muondataSim.Hits()->GetEntriesFast();
5707a463	220	AliMUONHit* mHit;
	221
	222	for (ihit=0; ihit<nhits; ihit++) {
b19dcb40	223	mHit = static_cast<AliMUONHit*>(muondataSim.Hits()->At(ihit));
5707a463	224	Int_t Nch = mHit->Chamber();
5707a463	225	Int_t hittrack = mHit->Track();
09f9a8f1	226	Float_t IdPart = mHit->Particle();
522b4724	227
5707a463	228	for (Int_t j=0;j<4;j++) {
5707a463	229	Int_t kch=11+j;
522b4724	230	if (hittrack==0) {
5707a463	231	if (Nch==kch && (IdPart==-13 \|\| IdPart==13) && NbHits[j]==0) NbHits[j]++;
522b4724	232	}
5707a463	233	}
5707a463	234	}
b19dcb40	235	muondataSim.ResetHits();
522b4724	236	} // end track loop
07d83718	237
5707a463	238	// 3/4 coincidence
5707a463	239	SumNbHits=NbHits[0]+NbHits[1]+NbHits[2]+NbHits[3];
522b4724	240
5707a463	241	if (SumNbHits==3 \|\| SumNbHits==4) {
	242	coincmuon++;
	243	ptcoinc->Fill(ptmu);
	244	}
	245
b19dcb40	246	} // end loop on event
	247
	248	if (coincmuon==0) {
	249	cout << " >>> <E> coincmuon = 0 after event loop " << "\n";
	250	cout << " >>> this probably means that input does not contain one (and only one) " << "\n";
	251	cout << " >>> muon track per event as it should " << "\n";
	252	cout << " >>> see README for further information " << "\n";
	253	cout << " >>> exiting now ! " << "\n";
	254	return;
	255	}
5707a463	256
b19dcb40	257
	258	MUONLoaderSim->UnloadHits();
	259	if (!readFromRP) {
	260	MUONLoaderSim->UnloadDigits();
	261	} else {
	262	MUONLoaderRec->UnloadRecPoints();
	263	}
	264	RunLoaderSim->UnloadKinematics();
5707a463	265
	266	fprintf(fdat,"\n");
	267	fprintf(fdat,"\n");
	268	fprintf(fdat," Number of events = %d \n",nevents);
	269	fprintf(fdat," Number of events with 3/4 coinc = %d \n",(Int_t)coincmuon);
5707a463	270	fprintf(fdat," Nomber of dimuons with 3/4 coinc Lpt cut = %d \n",(Int_t)lptmuon);
	271	fprintf(fdat," Number of dimuons with 3/4 coinc Hpt cut = %d \n",(Int_t)hptmuon);
	272	fprintf(fdat,"\n");
	273
522b4724	274	Double_t efficiency,error;
522b4724	275
5707a463	276	efficiency=lptmuon/coincmuon;
	277	error=efficiencyTMath::Sqrt((lptmuon+coincmuon)/(lptmuoncoincmuon));
	278	fprintf(fdat," Efficiency Lpt cut = %4.4f +/- %4.4f\n",efficiency,error);
	279
	280	efficiency=hptmuon/coincmuon;
	281	error=efficiencyTMath::Sqrt((hptmuon+coincmuon)/(hptmuoncoincmuon));
	282	fprintf(fdat," Efficiency Hpt cut = %4.4f +/- %4.4f\n",efficiency,error);
	283
	284	fclose(fdat);
	285
	286	Double_t x1,x2,xval,xerr;
	287
	288	for (Int_t i=0;i<50;i++) {
	289	x1=ptcoinc->GetBinContent(i+1);
	290
5707a463	291	x2=ptlpt->GetBinContent(i+1);
	292	if (x1!=0 && x2!=0) {
	293	xval=x2/x1;
	294	xerr=xvalTMath::Sqrt((x1+x2)/x1x2);
	295	ptlpt->SetBinContent(i+1,xval);
	296	ptlpt->SetBinError(i+1,0);
	297	} else {
	298	ptlpt->SetBinContent(i+1,0.);
	299	ptlpt->SetBinError(i+1,0.);
	300	}
	301
	302	x2=pthpt->GetBinContent(i+1);
	303	if (x1!=0 && x2!=0) {
	304	xval=x2/x1;
	305	xerr=xvalTMath::Sqrt((x1+x2)/x1x2);
	306	pthpt->SetBinContent(i+1,xval);
	307	pthpt->SetBinError(i+1,0);
	308	} else {
	309	pthpt->SetBinContent(i+1,0.);
	310	pthpt->SetBinError(i+1,0.);
	311	}
b19dcb40	312	}
5707a463	313
5707a463	314	TF1 *fitlpt = new TF1("fitlpt",fitArc,0.,5.,4);
	315	TF1 *fithpt = new TF1("fithpt",fitArc,0.,5.,4);
	316
	317	TCanvas *c1 = new TCanvas("c1","Trigger efficiency vs pt and y for single muon",200,0,900,400);
522b4724	318	c1->Divide(2,1);
5707a463	319
5707a463	320	c1->cd(1);
5707a463	321	ptlpt->SetMinimum(0.);
	322	ptlpt->SetMaximum(1.05);
	323	ptlpt->SetTitle("");
	324	ptlpt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
	325	ptlpt->GetYaxis()->SetTitle("Efficiency");
	326	//ptlpt->GetXaxis()->SetRange(3);
	327	ptlpt->Draw("");
	328	fitlpt->SetParameters(0.602,0.774,0.273,0.048);
	329	fitlpt->SetLineColor(2);
	330	fitlpt->SetLineWidth(2);
	331	fitlpt->Draw("SAME");
522b4724	332	TLegend * leg = new TLegend(0.5,0.38,0.85,0.53);
5707a463	333	leg = new TLegend(0.5,0.38,0.85,0.53);
	334	leg->SetBorderSize(0);
	335	leg->SetFillColor(0);
	336	leg->SetTextSize(0.05);
	337	leg->SetTextFont(22);
	338	leg->SetHeader("Lpt trigger pt cut");
	339	leg->AddEntry(fitlpt,"Ref.","l");
	340	leg->AddEntry(ptlpt,"New","l");
	341	leg->Draw("SAME");
	342
522b4724	343	c1->cd(2);
5707a463	344	pthpt->SetMinimum(0.);
	345	pthpt->SetMaximum(1.05);
	346	pthpt->SetTitle("");
	347	pthpt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
	348	pthpt->GetYaxis()->SetTitle("Efficiency");
	349	//pthpt->GetXaxis()->SetRange(3);
	350	pthpt->Draw("");
	351	fithpt->SetParameters(0.627,0.393,0.855,0.0081);
	352	fithpt->SetLineColor(2);
	353	fithpt->SetLineWidth(2);
	354	fithpt->Draw("SAME");
	355	leg = new TLegend(0.5,0.38,0.85,0.53);
	356	leg->SetBorderSize(0);
	357	leg->SetFillColor(0);
	358	leg->SetTextSize(0.05);
	359	leg->SetTextFont(22);
	360	leg->SetHeader("Hpt trigger pt cut");
	361	leg->AddEntry(fithpt,"Ref.","l");
	362	leg->AddEntry(pthpt,"New","l");
	363	leg->Draw("SAME");
	364
	365	c1->SaveAs("MUONTriggerEfficiencyPt.gif");
07d83718	366	c1->SaveAs("MUONTriggerEfficiencyPt.eps");
5707a463	367	}