[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 
   3 pt cuts. Results are compared to the reference (red curves).   
   To be used with (at least) 10000 events as follows
   AliGenBox * gener = new AliGenBox(1);
   gener->SetPtRange(0.,10.);
   gener->SetPhiRange(0., 360.);         
   gener->SetThetaRange(171.000,178.001);
   gener->SetPart(13);           // or -13
   gener->SetOrigin(0.,0., 0.);  
   gener->SetSigma(0.0, 0.0, 0.0);     

   Author: Fabien Guerin, LPC Clermont-Ferrand, Jan. 2006
*/


// 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"


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

// MUON includes
#include "AliMUON.h"
#include "AliMUONData.h"
#include "AliMUONHit.h"
#include "AliMUONChamber.h"
#include "AliMUONConstants.h"
#include "AliMUONDigit.h"
#include "AliMUONRawCluster.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 filename[10]="galice.root")
{

// 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 *ptapt    = new TH1F("ptapt","",50,0.15,5.); 
    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; 
    
    Double_t coincmuon=0;
    Double_t aptmuon=0;
    Double_t lptmuon=0;
    Double_t hptmuon=0;
    Double_t ptmu;

// output file
    char digitdat[100];
    sprintf(digitdat,"MUONTriggerEfficiencyPt.out");  
    FILE *fdat=fopen(digitdat,"w");

// Creating Run Loader and openning file containing Hits
    AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
    
    if (RunLoader ==0x0) {
	printf(">>> Error : Error Opening %s file \n",filename);
	return;
    }
    
    Int_t nevents = RunLoader->GetNumberOfEvents();
    
// loading hits
    AliLoader * MUONLoader = RunLoader->GetLoader("MUONLoader");     
    MUONLoader->LoadHits("READ");
    AliMUONData data_hits(MUONLoader,"MUON","MUON");
    
    TClonesArray * globalTrigger;
    AliMUONGlobalTrigger * gloTrg;

// loading trigger
    MUONLoader->LoadDigits("READ");  
    AliMUONData muondata(MUONLoader,"MUON","MUON");

// Loading kine
    RunLoader->LoadKinematics("READ");
    
    for (Int_t ievent=0; ievent<nevents; ievent++) {  // Event loop
	
	RunLoader->GetEvent(ievent);
	
	if (ievent%500==0) printf("ievent = %d \n",ievent);


// kine
	Int_t iparticle, nparticles;
	stack = RunLoader->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 
        muondata.SetTreeAddress("D,GLT");
        muondata.GetTriggerD();
    
        globalTrigger = muondata.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->SinglePlusApt()>=1 || gloTrg->SingleMinusApt()>=1 || gloTrg->SingleUndefApt()) {
		aptmuon++;
		ptapt->Fill(ptmu);
	    }
	    if (gloTrg->SinglePlusLpt()>=1 || gloTrg->SingleMinusLpt()>=1 || gloTrg->SingleUndefLpt()) {
		lptmuon++;
		ptlpt->Fill(ptmu);
	    }       
	    if (gloTrg->SinglePlusHpt()>=1 || gloTrg->SingleMinusHpt()>=1 || gloTrg->SingleUndefHpt()) {
		hptmuon++;
		pthpt->Fill(ptmu);
	    }
        } // end of loop on Global Trigger      
        muondata.ResetTrigger();  

// Hits
        RunLoader->GetEvent(ievent);  
        data_hits.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) data_hits.GetNtracks();      
     
        for (itrack=0; itrack<ntracks; itrack++) { // track loop
         data_hits.GetTrack(itrack); 

         Int_t ihit, nhits;
         nhits = (Int_t) data_hits.Hits()->GetEntriesFast();   
         AliMUONHit* mHit;
    
          for (ihit=0; ihit<nhits; ihit++) {
            mHit = static_cast<AliMUONHit*>(data_hits.Hits()->At(ihit));
            Int_t Nch        = mHit->Chamber(); 
            Int_t hittrack   = mHit->Track();
            Int_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]++; 
              }               
            }    
          }
          data_hits.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  
      
      MUONLoader->UnloadHits();
      MUONLoader->UnloadDigits();
      MUONLoader->UnloadRecPoints();
      RunLoader->UnloadKinematics();   

      delete RunLoader;
     
      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," Number of dimuons with 3/4 coinc Apt cut = %d \n",(Int_t)aptmuon);
      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=aptmuon/coincmuon;
      error=efficiency*TMath::Sqrt((aptmuon+coincmuon)/(aptmuon*coincmuon));
      fprintf(fdat," Efficiency Apt cut = %4.4f +/- %4.4f\n",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=ptapt->GetBinContent(i+1);
	  if (x1!=0 && x2!=0) {
	      xval=x2/x1;
	      xerr=xval*TMath::Sqrt((x1+x2)/x1*x2);   
	      ptapt->SetBinContent(i+1,xval);
	      ptapt->SetBinError(i+1,0);
	  } else {
	      ptapt->SetBinContent(i+1,0.);
	      ptapt->SetBinError(i+1,0.);     
	  }
	  
	  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 *fitapt = new TF1("fitapt",fitArch,0.,5.,4);     
      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(3,1);      
      
      c1->cd(1);
      ptapt->SetMinimum(0.);
      ptapt->SetMaximum(1.05);   
      ptapt->SetTitle("");
      ptapt->GetXaxis()->SetTitle("P_{T} (GeV/c)"); 
      ptapt->GetYaxis()->SetTitle("Efficiency");     
      //ptapt->GetXaxis()->SetRange(3);
      ptapt->Draw("");
      fitapt->SetParameters(1.853,1.57,-0.0203,-0.178);   
      fitapt->SetLineColor(2);
      fitapt->SetLineWidth(2);
      fitapt->Draw("SAME");    
      TLegend * 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("Apt trigger pt cut");
      leg->AddEntry(fitapt,"Ref.","l");
      leg->AddEntry(ptapt,"New","l");                 
      leg->Draw("SAME");     
      
      c1->cd(2);
      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");    
      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(3);
      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
	19	/* Macro to produce trigger single muon efficiency versus pt plots for the
	20	3 pt cuts. Results are compared to the reference (red curves).
	21	To be used with (at least) 10000 events as follows
	22	AliGenBox * gener = new AliGenBox(1);
	23	gener->SetPtRange(0.,10.);
	24	gener->SetPhiRange(0., 360.);
	25	gener->SetThetaRange(171.000,178.001);
	26	gener->SetPart(13); // or -13
	27	gener->SetOrigin(0.,0., 0.);
	28	gener->SetSigma(0.0, 0.0, 0.0);
	29
	30	Author: Fabien Guerin, LPC Clermont-Ferrand, Jan. 2006
	31	*/
	32
	33
	34	// ROOT includes
	35	#include "TBranch.h"
	36	#include "TClonesArray.h"
	37	#include "TFile.h"
	38	#include "TF1.h"
	39	#include "TH1.h"
	40	#include "TH2.h"
	41	#include "TParticle.h"
	42	#include "TTree.h"
	43	#include "TMath.h"
	44	#include "TStyle.h"
	45	#include "TCanvas.h"
	46	#include "TLegend.h"
	47
	48
	49	// STEER includes
	50	#include "AliRun.h"
	51	#include "AliDetector.h"
	52	#include "AliRunLoader.h"
	53	#include "AliHeader.h"
	54	#include "AliLoader.h"
	55	#include "AliStack.h"
	56	#include "AliSegmentation.h"
	57	#include "AliMC.h"
	58
	59	// MUON includes
	60	#include "AliMUON.h"
	61	#include "AliMUONData.h"
	62	#include "AliMUONHit.h"
	63	#include "AliMUONChamber.h"
	64	#include "AliMUONConstants.h"
65	#include "AliMUONDigit.h"
66	#include "AliMUONRawCluster.h"
67	#include "AliMUONGlobalTrigger.h"
68	#include "AliMUONLocalTrigger.h"
69	#include "AliMUONTrack.h"
70
71	Double_t fitArc(Double_t x,Double_t par)
72	{
73	Double_t h1=par[1]*(x[0]+par[2]);
74	return par[0]*TMath::ATan(TMath::Power(h1,6))+par[3];
75	}
76
77	Double_t fitArch(Double_t x,Double_t par)
78	{
79	Double_t h0=2*par[0]/(TMath::Pi());
80	Double_t h1=par[1]x[0]+par[2]x[0]*x[0];
81	return h0*TMath::TanH(h1)+par[3];
82	}
83
84	void MUONTriggerEfficiencyPt(char filename[10]="galice.root")
85	{
86
87	// define style
88	TStyle *st1 = new TStyle("st1","My STYLE");
89
90	st1->SetOptStat(0);
91	st1->SetOptFit(111);
92	st1->SetFrameFillColor(10);
93	st1->SetCanvasColor(10);
94	st1->SetFillColor(10);
95	st1->SetTitleBorderSize(0);
96	st1->SetTitleOffset(1.2,"XY");
97	st1->SetTitleSize(0.05,"XY");
98	st1->SetLabelSize(0.045,"XY");
99	st1->SetLabelFont(22,"XY");
100	st1->SetTitleFont(22,"XY");
101	st1->SetOptTitle(0);
102	st1->SetStatColor(10);
103	st1->SetStatFont(62);
104	st1->SetFitFormat("4.4g");
105	st1->SetPadLeftMargin(0.15);
106	st1->SetPadRightMargin(0.06);
107	st1->SetPadTopMargin(0.06);
108	st1->SetPadBottomMargin(0.15);
109	st1->cd();
110
111	gROOT->ForceStyle();
112	//TGaxis::SetMaxDigits(3);
113
114	// beginning of macro
115	TH1F *ptapt = new TH1F("ptapt","",50,0.15,5.);
116	TH1F *ptlpt = new TH1F("ptlpt","",50,0.15,5.);
117	TH1F *pthpt = new TH1F("pthpt","",50,0.15,5.);
118	TH1F *ptcoinc = new TH1F("ptcoinc","",50,0.15,5.);
119
120	TParticle *particle;
121	AliStack* stack;
122
123	Double_t coincmuon=0;
124	Double_t aptmuon=0;
125	Double_t lptmuon=0;
126	Double_t hptmuon=0;
127	Double_t ptmu;
128
129	// output file
130	char digitdat[100];
131	sprintf(digitdat,"MUONTriggerEfficiencyPt.out");
132	FILE *fdat=fopen(digitdat,"w");
133
134	// Creating Run Loader and openning file containing Hits
135	AliRunLoader * RunLoader = AliRunLoader::Open(filename,"MUONFolder","READ");
136
137	if (RunLoader ==0x0) {
138	printf(">>> Error : Error Opening %s file \n",filename);
139	return;
140	}
141
142	Int_t nevents = RunLoader->GetNumberOfEvents();
143
144	// loading hits
145	AliLoader * MUONLoader = RunLoader->GetLoader("MUONLoader");
146	MUONLoader->LoadHits("READ");
147	AliMUONData data_hits(MUONLoader,"MUON","MUON");
148
149	TClonesArray * globalTrigger;
150	AliMUONGlobalTrigger * gloTrg;
151
152	// loading trigger
153	MUONLoader->LoadDigits("READ");
154	AliMUONData muondata(MUONLoader,"MUON","MUON");
155
156	// Loading kine
157	RunLoader->LoadKinematics("READ");
158
159	for (Int_t ievent=0; ievent<nevents; ievent++) { // Event loop
160
161	RunLoader->GetEvent(ievent);
162
163	if (ievent%500==0) printf("ievent = %d \n",ievent);
164
165
166	// kine
167	Int_t iparticle, nparticles;
168	stack = RunLoader->Stack();
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
178	muondata.SetTreeAddress("D,GLT");
179	muondata.GetTriggerD();
180
181	globalTrigger = muondata.GlobalTrigger();
182
183	Int_t nglobals = (Int_t) globalTrigger->GetEntriesFast(); // should be 1
184	for (Int_t iglobal=0; iglobal<nglobals; iglobal++) { // Global Trigger
185	gloTrg = static_cast<AliMUONGlobalTrigger*>(globalTrigger->At(iglobal));
186
187	if (gloTrg->SinglePlusApt()>=1 \|\| gloTrg->SingleMinusApt()>=1 \|\| gloTrg->SingleUndefApt()) {
188	aptmuon++;
189	ptapt->Fill(ptmu);
190	}
191	if (gloTrg->SinglePlusLpt()>=1 \|\| gloTrg->SingleMinusLpt()>=1 \|\| gloTrg->SingleUndefLpt()) {
192	lptmuon++;
193	ptlpt->Fill(ptmu);
194	}
195	if (gloTrg->SinglePlusHpt()>=1 \|\| gloTrg->SingleMinusHpt()>=1 \|\| gloTrg->SingleUndefHpt()) {
196	hptmuon++;
197	pthpt->Fill(ptmu);
198	}
199	} // end of loop on Global Trigger
200	muondata.ResetTrigger();
201
202	// Hits
203	RunLoader->GetEvent(ievent);
204	data_hits.SetTreeAddress("H");
205
206	Int_t itrack, ntracks, NbHits[4];
207	Int_t SumNbHits;
208
209	for (Int_t j=0; j<4; j++) NbHits[j]=0;
210
211	ntracks = (Int_t) data_hits.GetNtracks();
212
213	for (itrack=0; itrack<ntracks; itrack++) { // track loop
214	data_hits.GetTrack(itrack);
215
216	Int_t ihit, nhits;
217	nhits = (Int_t) data_hits.Hits()->GetEntriesFast();
218	AliMUONHit* mHit;
219
220	for (ihit=0; ihit<nhits; ihit++) {
221	mHit = static_cast<AliMUONHit*>(data_hits.Hits()->At(ihit));
222	Int_t Nch = mHit->Chamber();
223	Int_t hittrack = mHit->Track();
224	Int_t IdPart = mHit->Particle();
225
226	for (Int_t j=0;j<4;j++) {
227	Int_t kch=11+j;
228	if (hittrack==0) {
229	if (Nch==kch && (IdPart==-13 \|\| IdPart==13) && NbHits[j]==0) NbHits[j]++;
230	}
231	}
232	}
233	data_hits.ResetHits();
234	} // end track loop
235
236
237	// 3/4 coincidence
238	SumNbHits=NbHits[0]+NbHits[1]+NbHits[2]+NbHits[3];
239
240	if (SumNbHits==3 \|\| SumNbHits==4) {
241	coincmuon++;
242	ptcoinc->Fill(ptmu);
243	}
244
245	} // end loop on event
246
247	MUONLoader->UnloadHits();
248	MUONLoader->UnloadDigits();
249	MUONLoader->UnloadRecPoints();
250	RunLoader->UnloadKinematics();
251
252	delete RunLoader;
253
254	fprintf(fdat,"\n");
255	fprintf(fdat,"\n");
256	fprintf(fdat," Number of events = %d \n",nevents);
257	fprintf(fdat," Number of events with 3/4 coinc = %d \n",(Int_t)coincmuon);
258	fprintf(fdat," Number of dimuons with 3/4 coinc Apt cut = %d \n",(Int_t)aptmuon);
259	fprintf(fdat," Nomber of dimuons with 3/4 coinc Lpt cut = %d \n",(Int_t)lptmuon);
260	fprintf(fdat," Number of dimuons with 3/4 coinc Hpt cut = %d \n",(Int_t)hptmuon);
261	fprintf(fdat,"\n");
262
263	Double_t efficiency,error;
264
265	efficiency=aptmuon/coincmuon;
266	error=efficiencyTMath::Sqrt((aptmuon+coincmuon)/(aptmuoncoincmuon));
267	fprintf(fdat," Efficiency Apt cut = %4.4f +/- %4.4f\n",efficiency,error);
268
269	efficiency=lptmuon/coincmuon;
270	error=efficiencyTMath::Sqrt((lptmuon+coincmuon)/(lptmuoncoincmuon));
271	fprintf(fdat," Efficiency Lpt cut = %4.4f +/- %4.4f\n",efficiency,error);
272
273	efficiency=hptmuon/coincmuon;
274	error=efficiencyTMath::Sqrt((hptmuon+coincmuon)/(hptmuoncoincmuon));
275	fprintf(fdat," Efficiency Hpt cut = %4.4f +/- %4.4f\n",efficiency,error);
276
277	fclose(fdat);
278
279	Double_t x1,x2,xval,xerr;
280
281	for (Int_t i=0;i<50;i++) {
282	x1=ptcoinc->GetBinContent(i+1);
283
284	x2=ptapt->GetBinContent(i+1);
285	if (x1!=0 && x2!=0) {
286	xval=x2/x1;
287	xerr=xvalTMath::Sqrt((x1+x2)/x1x2);
288	ptapt->SetBinContent(i+1,xval);
289	ptapt->SetBinError(i+1,0);
290	} else {
291	ptapt->SetBinContent(i+1,0.);
292	ptapt->SetBinError(i+1,0.);
293	}
294
295	x2=ptlpt->GetBinContent(i+1);
296	if (x1!=0 && x2!=0) {
297	xval=x2/x1;
298	xerr=xvalTMath::Sqrt((x1+x2)/x1x2);
299	ptlpt->SetBinContent(i+1,xval);
300	ptlpt->SetBinError(i+1,0);
301	} else {
302	ptlpt->SetBinContent(i+1,0.);
303	ptlpt->SetBinError(i+1,0.);
304	}
305
306	x2=pthpt->GetBinContent(i+1);
307	if (x1!=0 && x2!=0) {
308	xval=x2/x1;
309	xerr=xvalTMath::Sqrt((x1+x2)/x1x2);
310	pthpt->SetBinContent(i+1,xval);
311	pthpt->SetBinError(i+1,0);
312	} else {
313	pthpt->SetBinContent(i+1,0.);
314	pthpt->SetBinError(i+1,0.);
315	}
316	}
317
318
319	TF1 *fitapt = new TF1("fitapt",fitArch,0.,5.,4);
320	TF1 *fitlpt = new TF1("fitlpt",fitArc,0.,5.,4);
321	TF1 *fithpt = new TF1("fithpt",fitArc,0.,5.,4);
322
323	TCanvas *c1 = new TCanvas("c1","Trigger efficiency vs pt and y for single muon",200,0,900,400);
324	c1->Divide(3,1);
325
326	c1->cd(1);
327	ptapt->SetMinimum(0.);
328	ptapt->SetMaximum(1.05);
329	ptapt->SetTitle("");
330	ptapt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
331	ptapt->GetYaxis()->SetTitle("Efficiency");
332	//ptapt->GetXaxis()->SetRange(3);
333	ptapt->Draw("");
334	fitapt->SetParameters(1.853,1.57,-0.0203,-0.178);
335	fitapt->SetLineColor(2);
336	fitapt->SetLineWidth(2);
337	fitapt->Draw("SAME");
338	TLegend * leg = new TLegend(0.5,0.38,0.85,0.53);
339	leg->SetBorderSize(0);
340	leg->SetFillColor(0);
341	leg->SetTextSize(0.05);
342	leg->SetTextFont(22);
343	leg->SetHeader("Apt trigger pt cut");
344	leg->AddEntry(fitapt,"Ref.","l");
345	leg->AddEntry(ptapt,"New","l");
346	leg->Draw("SAME");
347
348	c1->cd(2);
349	ptlpt->SetMinimum(0.);
350	ptlpt->SetMaximum(1.05);
351	ptlpt->SetTitle("");
352	ptlpt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
353	ptlpt->GetYaxis()->SetTitle("Efficiency");
354	//ptlpt->GetXaxis()->SetRange(3);
355	ptlpt->Draw("");
356	fitlpt->SetParameters(0.602,0.774,0.273,0.048);
357	fitlpt->SetLineColor(2);
358	fitlpt->SetLineWidth(2);
359	fitlpt->Draw("SAME");
360	leg = new TLegend(0.5,0.38,0.85,0.53);
361	leg->SetBorderSize(0);
362	leg->SetFillColor(0);
363	leg->SetTextSize(0.05);
364	leg->SetTextFont(22);
365	leg->SetHeader("Lpt trigger pt cut");
366	leg->AddEntry(fitlpt,"Ref.","l");
367	leg->AddEntry(ptlpt,"New","l");
368	leg->Draw("SAME");
369
370	c1->cd(3);
371	pthpt->SetMinimum(0.);
372	pthpt->SetMaximum(1.05);
373	pthpt->SetTitle("");
374	pthpt->GetXaxis()->SetTitle("P_{T} (GeV/c)");
375	pthpt->GetYaxis()->SetTitle("Efficiency");
376	//pthpt->GetXaxis()->SetRange(3);
377	pthpt->Draw("");
378	fithpt->SetParameters(0.627,0.393,0.855,0.0081);
379	fithpt->SetLineColor(2);
380	fithpt->SetLineWidth(2);
381	fithpt->Draw("SAME");
382	leg = new TLegend(0.5,0.38,0.85,0.53);
383	leg->SetBorderSize(0);
384	leg->SetFillColor(0);
385	leg->SetTextSize(0.05);
386	leg->SetTextFont(22);
387	leg->SetHeader("Hpt trigger pt cut");
388	leg->AddEntry(fithpt,"Ref.","l");
389	leg->AddEntry(pthpt,"New","l");
390	leg->Draw("SAME");
391
392	c1->SaveAs("MUONTriggerEfficiencyPt.gif");
393	c1->SaveAs("MUONTriggerEfficiencyPt.eps");
394
395	}