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

// macro to make invariant mass plots
// for combinations of 2 muons with opposite charges,
// from root file "MUONtrackReco.root" containing the result of track reconstruction,
// generated by the macro "MUONrecoNtuple.C".
// Histograms are stored on the "MUONmassPlot.root" file.
// A model for macros using the Ntuple in the file "MUONtrackReco.root"
// may be found in "MUONtrackRecoModel.C":
// it has been obtained by reading with Root a file "MUONtrackReco.root",
// and executing the command:
// MUONtrackReco->MakeCode("MUONtrackRecoModel.C")

// Arguments:
//   FirstEvent (default 0)
//   LastEvent (default 0)
//   ResType (default 553)
//      553 for Upsilon, anything else for J/Psi
//   Chi2Cut (default 100)
//      to keep only tracks with chi2 per d.o.f. < Chi2Cut
//   PtCutMin (default 1)
//      to keep only tracks with transverse momentum > PtCutMin
//   PtCutMax (default 10000)
//      to keep only tracks with transverse momentum < PtCutMax
//   massMin (default 9.17 for Upsilon) 
//      &  massMax (default 9.77 for Upsilon) 
//         to calculate the reconstruction efficiency for resonances with invariant mass
//         massMin < mass < massMax.

// compile MUONrecoNtuple.C and load shared library in this macro
// Add parameters and histograms for analysis 

void MUONmassPlot(Int_t FirstEvent = 0, Int_t LastEvent = 0, Int_t ResType = 553, Float_t Chi2Cut = 100., Float_t PtCutMin = 1., Float_t PtCutMax = 10000.,Float_t massMin = 9.17,Float_t massMax = 9.77)
{
  cout << "MUONmassPlot " << endl;
  cout << "FirstEvent " << FirstEvent << endl;
  cout << "LastEvent " << LastEvent << endl;
  cout << "ResType " << ResType << endl;
//   cout << "Nsig " << Nsig << endl;
  cout << "Chi2Cut " << Chi2Cut << endl;
  cout << "PtCutMin " << PtCutMin << endl;
  cout << "PtCutMax " << PtCutMax << endl;
  cout << "massMin " << massMin << endl;
  cout << "massMax " << massMax << endl;


  //////////////////////////////////////////////////////////
  //   This file has been automatically generated 
  //     (Thu Sep 21 14:53:11 2000 by ROOT version2.25/02)
  //   from TTree MUONtrackReco/MUONtrackReco
  //   found on file: MUONtrackReco.root
  //////////////////////////////////////////////////////////


  //Reset ROOT and connect tree file
  gROOT->Reset();
// Dynamically link some shared libs

  if (gClassTable->GetID("AliRun") < 0) {
    gSystem->SetIncludePath("-I$ALICE_ROOT/MUON -I$ALICE_ROOT/STEER");
    gROOT->LoadMacro("loadlibs.C");
    loadlibs();
// compile MUONrecoNtuple and load shared library 
    gSystem->CompileMacro("$(ALICE_ROOT)/macros/MUONrecoNtuple.C","kf");
  }

  TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("MUONtrackReco.root");
  if (!f) {
    f = new TFile("MUONtrackReco.root");
  }
  TTree *MUONtrackReco = (TTree*)gDirectory->Get("MUONtrackReco");
  MUONtrackReco->SetMakeClass(1);                                      

//Declaration of leaves types
  Int_t           fEvent;
  UInt_t          fUniqueID;
  UInt_t          fBits;
  Int_t           Tracks_;
  Int_t           Tracks_fCharge[500];
  Float_t         Tracks_fPxRec[500];
  Float_t         Tracks_fPyRec[500];
  Float_t         Tracks_fPzRec[500];
  Float_t         Tracks_fZRec[500];
  Float_t         Tracks_fZRec1[500];
  Int_t           Tracks_fNHits[500];
  Float_t         Tracks_fChi2[500];
  Float_t         Tracks_fPxGen[500];
  Float_t         Tracks_fPyGen[500];
  Float_t         Tracks_fPzGen[500];
  UInt_t          Tracks_fUniqueID[500];
  UInt_t          Tracks_fBits[500];

  //Set branch addresses
  //MUONtrackReco->SetBranchAddress("Header",&Header);
  MUONtrackReco->SetBranchAddress("fEvent",&fEvent);
  MUONtrackReco->SetBranchAddress("fUniqueID",&fUniqueID);
  MUONtrackReco->SetBranchAddress("fBits",&fBits);
//   MUONtrackReco->SetBranchAddress("Tracks_",&Tracks_);
  MUONtrackReco->SetBranchAddress("Tracks",&Tracks_);
  MUONtrackReco->SetBranchAddress("Tracks.fCharge",Tracks_fCharge);
  MUONtrackReco->SetBranchAddress("Tracks.fPxRec",Tracks_fPxRec);
  MUONtrackReco->SetBranchAddress("Tracks.fPyRec",Tracks_fPyRec);
  MUONtrackReco->SetBranchAddress("Tracks.fPzRec",Tracks_fPzRec);
  MUONtrackReco->SetBranchAddress("Tracks.fZRec",Tracks_fZRec);
  MUONtrackReco->SetBranchAddress("Tracks.fZRec1",Tracks_fZRec1);
  MUONtrackReco->SetBranchAddress("Tracks.fNHits",Tracks_fNHits);
  MUONtrackReco->SetBranchAddress("Tracks.fChi2",Tracks_fChi2);
  MUONtrackReco->SetBranchAddress("Tracks.fPxGen",Tracks_fPxGen);
  MUONtrackReco->SetBranchAddress("Tracks.fPyGen",Tracks_fPyGen);
  MUONtrackReco->SetBranchAddress("Tracks.fPzGen",Tracks_fPzGen);
  MUONtrackReco->SetBranchAddress("Tracks.fUniqueID",Tracks_fUniqueID);
  MUONtrackReco->SetBranchAddress("Tracks.fBits",Tracks_fBits);

//     This is the loop skeleton
//       To read only selected branches, Insert statements like:
// MUONtrackReco->SetBranchStatus("*",0);  // disable all branches
// TTreePlayer->SetBranchStatus("branchname",1);  // activate branchname

  Int_t nentries = MUONtrackReco->GetEntries();

  Int_t nbytes = 0;
  //   for (Int_t i=0; i<nentries;i++) {
  //      nbytes += MUONtrackReco->GetEntry(i);
  //   }

  /////////////////////////////////////////////////////////////////
  // Here comes the specialized part for MUONmassPlot
  /////////////////////////////////////////////////////////////////

  // File for histograms and histogram booking
  TFile *histoFile = new TFile("MUONmassPlot.root", "RECREATE");
  TH1F *hPtMuon = new TH1F("hPtMuon", "Muon Pt (GeV/c)", 100, 0., 20.);
  TH1F *hPMuon = new TH1F("hPMuon", "Muon P (GeV/c)", 100, 0., 200.);
  TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "Muon track chi2/d.o.f.", 100, 0., 20.);
  TH1F *hInvMassAll = new TH1F("hInvMassAll", "Mu+Mu- invariant mass (GeV/c2)", 480, 0., 12.);
  if (ResType == 553) {
    TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around Upsilon", 60, 8., 11.);
  } else {
    TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around J/Psi", 80, 0., 5.);
  }

  TH1F *hNumberOfTrack = new TH1F("hNumberOfTrack","nb of track /evt ",20,-0.5,19.5);
  TH1F *hRapMuon = new TH1F("hRapMuon"," Muon Rapidity",50,2.,4.5);
  TH1F *hRapResonance = new TH1F("hRapResonance"," Resonance Rapidity",50,2.,4.5);
  TH1F *hPtResonance = new TH1F("hPtResonance", "Resonance Pt (GeV/c)", 100, 0., 20.);
  Int_t EventInMass = 0;
  Float_t muonMass = 0.105658389;
  Float_t UpsilonMass = 9.46037;
  Float_t JPsiMass = 3.097;
  // Loop over events
  for (Int_t event = FirstEvent; event <= TMath::Min(LastEvent, nentries - 1); event++) {
    // get current event
    Int_t NumberOfTrack = 0;
    cout <<" event: " << event <<endl;
    cout << " number of tracks: " << Tracks_ <<endl;
    nbytes += MUONtrackReco->GetEntry(event);
    // loop over all reconstructed tracks (also first track of combination)
    for (Int_t t1 = 0; t1 < Tracks_; t1++) {
      // transverse momentum
      Float_t pt1 = TMath::Sqrt(Tracks_fPxRec[t1] * Tracks_fPxRec[t1] + Tracks_fPyRec[t1] * Tracks_fPyRec[t1]);

      // total momentum
      Float_t p1 = TMath::Sqrt(Tracks_fPxRec[t1] * Tracks_fPxRec[t1] + Tracks_fPyRec[t1] * Tracks_fPyRec[t1] + Tracks_fPzRec[t1] * Tracks_fPzRec[t1] );
      // Rapidity
      Float_t rapMuon1 = rapParticle(Tracks_fPxRec[t1],Tracks_fPyRec[t1],Tracks_fPzRec[t1],muonMass);

      // chi2 per d.o.f.
      Float_t ch1 = Tracks_fChi2[t1] / (2.0 * Tracks_fNHits[t1] - 5);
      printf(" px %f py %f pz %f NHits %d  Norm.chi2 %f \n",Tracks_fPxRec[t1],Tracks_fPyRec[t1],Tracks_fPzRec[t1],Tracks_fNHits[t1],ch1);
      // condition for good track (Chi2Cut and PtCut)
      if ((ch1 < Chi2Cut) && (pt1 > PtCutMin) && (pt1 < PtCutMax)) {
	// fill histos hPtMuon and hChi2PerDof
	NumberOfTrack++;
	hPtMuon->Fill(pt1);
	hPMuon->Fill(p1);
	hChi2PerDof->Fill(ch1);
	hRapMuon->Fill(rapMuon1);
	// loop over second track of combination
	for (Int_t t2 = t1 + 1; t2 < Tracks_; t2++) {
	  // transverse momentum
	  Float_t pt2 = TMath::Sqrt(Tracks_fPxRec[t2] * Tracks_fPxRec[t2] + Tracks_fPyRec[t2] * Tracks_fPyRec[t2]);
	  // chi2 per d.o.f.
	  Float_t ch2 = Tracks_fChi2[t2] / (2.0 * Tracks_fNHits[t2] - 5);
	  // condition for good track (Chi2Cut and PtCut)
	  if ((ch2 < Chi2Cut) && (pt2 > PtCutMin)  && (pt2 < PtCutMax)) {
	    // condition for opposite charges
	    if ((Tracks_fCharge[t1] * Tracks_fCharge[t2]) == -1) {
	      // invariant mass
	      Float_t invMass = MuPlusMuMinusMass(Tracks_fPxRec[t1], Tracks_fPyRec[t1], Tracks_fPzRec[t1], Tracks_fPxRec[t2], Tracks_fPyRec[t2], Tracks_fPzRec[t2],muonMass);
	      // fill histos hInvMassAll and hInvMassRes
	      hInvMassAll->Fill(invMass);
	      hInvMassRes->Fill(invMass);
	      Float_t pxRes = Tracks_fPxRec[t1]+Tracks_fPxRec[t2];
	      Float_t pyRes = Tracks_fPyRec[t1]+Tracks_fPyRec[t2];
	      Float_t pzRes = Tracks_fPzRec[t1]+Tracks_fPzRec[t2];
	      if (invMass > massMin && invMass < massMax) {
		EventInMass++;
		Float_t rapRes  = 0;
		if (ResType == 553) { 
		  rapRes  =  rapParticle(pxRes,pyRes,pzRes,UpsilonMass); 
		} else {
		  rapRes  =  rapParticle(pxRes,pyRes,pzRes,JPsiMass); 
		}
  		hRapResonance->Fill(rapRes);
  		hPtResonance->Fill(TMath::Sqrt(pxRes*pxRes+pyRes*pyRes));
	      }

	    } //if ((Tracks_fCharge[t1] * Tracks_fCharge[t2]) == -1)
	  } // if ((Tracks_fChi2[t2] < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax))
	} // for (Int_t t2 = t1 + 1; t2 < Tracks_; t2++)
      } // if ((Tracks_fChi2[t1] < Chi2Cut) && (pt1 > PtCutMin)&& (pt1 < PtCutMax) )
    } // for (Int_t t1 = 0; t1 < Tracks_; t1++)
    hNumberOfTrack->Fill(NumberOfTrack);
  } // for (Int_t event = FirstEvent;

  histoFile->Write();
  histoFile->Close();

  cout << "MUONmassPlot " << endl;
  cout << "FirstEvent " << FirstEvent << endl;
  cout << "LastEvent " << LastEvent << endl;
  cout << "ResType " << ResType << endl;
//   cout << "Nsig " << Nsig << endl;
  cout << "Chi2Cut " << Chi2Cut << endl;
  cout << "PtCutMin " << PtCutMin << endl;
  cout << "PtCutMax " << PtCutMax << endl;
  cout << "massMin " << massMin << endl;
  cout << "massMax " << massMax << endl;
  cout << "EventInMass " << EventInMass << endl;
}


Float_t MuPlusMuMinusMass(Float_t Px1, Float_t Py1, Float_t Pz1, Float_t Px2, Float_t Py2, Float_t Pz2, Float_t muonMass)
{
  // return invariant mass for particle 1 & 2 whose masses are equal to muonMass
 
  Float_t e1 = TMath::Sqrt(muonMass * muonMass + Px1 * Px1 + Py1 * Py1 + Pz1 * Pz1);
  Float_t e2 = TMath::Sqrt(muonMass * muonMass + Px2 * Px2 + Py2 * Py2 + Pz2 * Pz2);
  return (TMath::Sqrt(2.0 * (muonMass * muonMass + e1 * e2 - Px1 * Px2 - Py1 * Py2 - Pz1 * Pz2)));
}

Float_t rapParticle(Float_t Px, Float_t Py, Float_t Pz, Float_t Mass)
{
  // return rapidity for particle 

  // Particle energy
  Float_t Energy = TMath::Sqrt(Px*Px+Py*Py+Pz*Pz+Mass*Mass);
  // Rapidity
  Float_t Rapidity = 0.5*TMath::Log((Energy+Pz)/(Energy-Pz));
  return Rapidity;
}
Commit	Line	Data
c2a319d4	1	// macro to make invariant mass plots
	2	// for combinations of 2 muons with opposite charges,
	3	// from root file "MUONtrackReco.root" containing the result of track reconstruction,
	4	// generated by the macro "MUONrecoNtuple.C".
	5	// Histograms are stored on the "MUONmassPlot.root" file.
	6	// A model for macros using the Ntuple in the file "MUONtrackReco.root"
	7	// may be found in "MUONtrackRecoModel.C":
	8	// it has been obtained by reading with Root a file "MUONtrackReco.root",
	9	// and executing the command:
	10	// MUONtrackReco->MakeCode("MUONtrackRecoModel.C")
	11
	12	// Arguments:
	13	// FirstEvent (default 0)
	14	// LastEvent (default 0)
	15	// ResType (default 553)
	16	// 553 for Upsilon, anything else for J/Psi
c2a319d4	17	// Chi2Cut (default 100)
c2a319d4	18	// to keep only tracks with chi2 per d.o.f. < Chi2Cut
8904b7e3	19	// PtCutMin (default 1)
	20	// to keep only tracks with transverse momentum > PtCutMin
	21	// PtCutMax (default 10000)
	22	// to keep only tracks with transverse momentum < PtCutMax
	23	// massMin (default 9.17 for Upsilon)
	24	// & massMax (default 9.77 for Upsilon)
	25	// to calculate the reconstruction efficiency for resonances with invariant mass
	26	// massMin < mass < massMax.
	27
	28	// compile MUONrecoNtuple.C and load shared library in this macro
	29	// Add parameters and histograms for analysis
	30
	31	void MUONmassPlot(Int_t FirstEvent = 0, Int_t LastEvent = 0, Int_t ResType = 553, Float_t Chi2Cut = 100., Float_t PtCutMin = 1., Float_t PtCutMax = 10000.,Float_t massMin = 9.17,Float_t massMax = 9.77)
c2a319d4	32	{
8904b7e3	33	cout << "MUONmassPlot " << endl;
	34	cout << "FirstEvent " << FirstEvent << endl;
	35	cout << "LastEvent " << LastEvent << endl;
	36	cout << "ResType " << ResType << endl;
	37	// cout << "Nsig " << Nsig << endl;
	38	cout << "Chi2Cut " << Chi2Cut << endl;
	39	cout << "PtCutMin " << PtCutMin << endl;
	40	cout << "PtCutMax " << PtCutMax << endl;
	41	cout << "massMin " << massMin << endl;
	42	cout << "massMax " << massMax << endl;
	43
c2a319d4	44
	45	//////////////////////////////////////////////////////////
	46	// This file has been automatically generated
	47	// (Thu Sep 21 14:53:11 2000 by ROOT version2.25/02)
	48	// from TTree MUONtrackReco/MUONtrackReco
	49	// found on file: MUONtrackReco.root
	50	//////////////////////////////////////////////////////////
	51
	52
	53	//Reset ROOT and connect tree file
	54	gROOT->Reset();
8904b7e3	55	// Dynamically link some shared libs
	56
	57	if (gClassTable->GetID("AliRun") < 0) {
	58	gSystem->SetIncludePath("-I$ALICE_ROOT/MUON -I$ALICE_ROOT/STEER");
	59	gROOT->LoadMacro("loadlibs.C");
	60	loadlibs();
	61	// compile MUONrecoNtuple and load shared library
	62	gSystem->CompileMacro("$(ALICE_ROOT)/macros/MUONrecoNtuple.C","kf");
	63	}
	64
c2a319d4	65	TFile f = (TFile)gROOT->GetListOfFiles()->FindObject("MUONtrackReco.root");
	66	if (!f) {
	67	f = new TFile("MUONtrackReco.root");
	68	}
	69	TTree MUONtrackReco = (TTree)gDirectory->Get("MUONtrackReco");
eff57104	70	MUONtrackReco->SetMakeClass(1);
c2a319d4	71
	72	//Declaration of leaves types
	73	Int_t fEvent;
	74	UInt_t fUniqueID;
	75	UInt_t fBits;
	76	Int_t Tracks_;
8904b7e3	77	Int_t Tracks_fCharge[500];
	78	Float_t Tracks_fPxRec[500];
	79	Float_t Tracks_fPyRec[500];
	80	Float_t Tracks_fPzRec[500];
	81	Float_t Tracks_fZRec[500];
	82	Float_t Tracks_fZRec1[500];
	83	Int_t Tracks_fNHits[500];
	84	Float_t Tracks_fChi2[500];
	85	Float_t Tracks_fPxGen[500];
	86	Float_t Tracks_fPyGen[500];
	87	Float_t Tracks_fPzGen[500];
	88	UInt_t Tracks_fUniqueID[500];
	89	UInt_t Tracks_fBits[500];
c2a319d4	90
	91	//Set branch addresses
	92	//MUONtrackReco->SetBranchAddress("Header",&Header);
	93	MUONtrackReco->SetBranchAddress("fEvent",&fEvent);
	94	MUONtrackReco->SetBranchAddress("fUniqueID",&fUniqueID);
	95	MUONtrackReco->SetBranchAddress("fBits",&fBits);
eff57104	96	// MUONtrackReco->SetBranchAddress("Tracks_",&Tracks_);
eff57104	97	MUONtrackReco->SetBranchAddress("Tracks",&Tracks_);
c2a319d4	98	MUONtrackReco->SetBranchAddress("Tracks.fCharge",Tracks_fCharge);
	99	MUONtrackReco->SetBranchAddress("Tracks.fPxRec",Tracks_fPxRec);
	100	MUONtrackReco->SetBranchAddress("Tracks.fPyRec",Tracks_fPyRec);
	101	MUONtrackReco->SetBranchAddress("Tracks.fPzRec",Tracks_fPzRec);
	102	MUONtrackReco->SetBranchAddress("Tracks.fZRec",Tracks_fZRec);
	103	MUONtrackReco->SetBranchAddress("Tracks.fZRec1",Tracks_fZRec1);
	104	MUONtrackReco->SetBranchAddress("Tracks.fNHits",Tracks_fNHits);
	105	MUONtrackReco->SetBranchAddress("Tracks.fChi2",Tracks_fChi2);
	106	MUONtrackReco->SetBranchAddress("Tracks.fPxGen",Tracks_fPxGen);
	107	MUONtrackReco->SetBranchAddress("Tracks.fPyGen",Tracks_fPyGen);
	108	MUONtrackReco->SetBranchAddress("Tracks.fPzGen",Tracks_fPzGen);
	109	MUONtrackReco->SetBranchAddress("Tracks.fUniqueID",Tracks_fUniqueID);
	110	MUONtrackReco->SetBranchAddress("Tracks.fBits",Tracks_fBits);
	111
	112	// This is the loop skeleton
	113	// To read only selected branches, Insert statements like:
	114	// MUONtrackReco->SetBranchStatus("*",0); // disable all branches
	115	// TTreePlayer->SetBranchStatus("branchname",1); // activate branchname
	116
	117	Int_t nentries = MUONtrackReco->GetEntries();
	118
	119	Int_t nbytes = 0;
	120	// for (Int_t i=0; i<nentries;i++) {
	121	// nbytes += MUONtrackReco->GetEntry(i);
	122	// }
	123
	124	/////////////////////////////////////////////////////////////////
	125	// Here comes the specialized part for MUONmassPlot
	126	/////////////////////////////////////////////////////////////////
	127
	128	// File for histograms and histogram booking
	129	TFile *histoFile = new TFile("MUONmassPlot.root", "RECREATE");
	130	TH1F *hPtMuon = new TH1F("hPtMuon", "Muon Pt (GeV/c)", 100, 0., 20.);
8904b7e3	131	TH1F *hPMuon = new TH1F("hPMuon", "Muon P (GeV/c)", 100, 0., 200.);
c2a319d4	132	TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "Muon track chi2/d.o.f.", 100, 0., 20.);
8904b7e3	133	TH1F *hInvMassAll = new TH1F("hInvMassAll", "Mu+Mu- invariant mass (GeV/c2)", 480, 0., 12.);
	134	if (ResType == 553) {
	135	TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around Upsilon", 60, 8., 11.);
	136	} else {
	137	TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around J/Psi", 80, 0., 5.);
	138	}
c2a319d4	139
8904b7e3	140	TH1F *hNumberOfTrack = new TH1F("hNumberOfTrack","nb of track /evt ",20,-0.5,19.5);
	141	TH1F *hRapMuon = new TH1F("hRapMuon"," Muon Rapidity",50,2.,4.5);
	142	TH1F *hRapResonance = new TH1F("hRapResonance"," Resonance Rapidity",50,2.,4.5);
	143	TH1F *hPtResonance = new TH1F("hPtResonance", "Resonance Pt (GeV/c)", 100, 0., 20.);
	144	Int_t EventInMass = 0;
	145	Float_t muonMass = 0.105658389;
	146	Float_t UpsilonMass = 9.46037;
	147	Float_t JPsiMass = 3.097;
c2a319d4	148	// Loop over events
	149	for (Int_t event = FirstEvent; event <= TMath::Min(LastEvent, nentries - 1); event++) {
	150	// get current event
8904b7e3	151	Int_t NumberOfTrack = 0;
	152	cout <<" event: " << event <<endl;
	153	cout << " number of tracks: " << Tracks_ <<endl;
c2a319d4	154	nbytes += MUONtrackReco->GetEntry(event);
	155	// loop over all reconstructed tracks (also first track of combination)
	156	for (Int_t t1 = 0; t1 < Tracks_; t1++) {
	157	// transverse momentum
	158	Float_t pt1 = TMath::Sqrt(Tracks_fPxRec[t1] * Tracks_fPxRec[t1] + Tracks_fPyRec[t1] * Tracks_fPyRec[t1]);
8904b7e3	159
	160	// total momentum
	161	Float_t p1 = TMath::Sqrt(Tracks_fPxRec[t1] * Tracks_fPxRec[t1] + Tracks_fPyRec[t1] * Tracks_fPyRec[t1] + Tracks_fPzRec[t1] * Tracks_fPzRec[t1] );
	162	// Rapidity
	163	Float_t rapMuon1 = rapParticle(Tracks_fPxRec[t1],Tracks_fPyRec[t1],Tracks_fPzRec[t1],muonMass);
	164
c2a319d4	165	// chi2 per d.o.f.
c2a319d4	166	Float_t ch1 = Tracks_fChi2[t1] / (2.0 * Tracks_fNHits[t1] - 5);
8904b7e3	167	printf(" px %f py %f pz %f NHits %d Norm.chi2 %f \n",Tracks_fPxRec[t1],Tracks_fPyRec[t1],Tracks_fPzRec[t1],Tracks_fNHits[t1],ch1);
c2a319d4	168	// condition for good track (Chi2Cut and PtCut)
8904b7e3	169	if ((ch1 < Chi2Cut) && (pt1 > PtCutMin) && (pt1 < PtCutMax)) {
c2a319d4	170	// fill histos hPtMuon and hChi2PerDof
8904b7e3	171	NumberOfTrack++;
c2a319d4	172	hPtMuon->Fill(pt1);
8904b7e3	173	hPMuon->Fill(p1);
c2a319d4	174	hChi2PerDof->Fill(ch1);
8904b7e3	175	hRapMuon->Fill(rapMuon1);
c2a319d4	176	// loop over second track of combination
	177	for (Int_t t2 = t1 + 1; t2 < Tracks_; t2++) {
	178	// transverse momentum
	179	Float_t pt2 = TMath::Sqrt(Tracks_fPxRec[t2] * Tracks_fPxRec[t2] + Tracks_fPyRec[t2] * Tracks_fPyRec[t2]);
	180	// chi2 per d.o.f.
	181	Float_t ch2 = Tracks_fChi2[t2] / (2.0 * Tracks_fNHits[t2] - 5);
	182	// condition for good track (Chi2Cut and PtCut)
8904b7e3	183	if ((ch2 < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax)) {
c2a319d4	184	// condition for opposite charges
	185	if ((Tracks_fCharge[t1] * Tracks_fCharge[t2]) == -1) {
	186	// invariant mass
8904b7e3	187	Float_t invMass = MuPlusMuMinusMass(Tracks_fPxRec[t1], Tracks_fPyRec[t1], Tracks_fPzRec[t1], Tracks_fPxRec[t2], Tracks_fPyRec[t2], Tracks_fPzRec[t2],muonMass);
c2a319d4	188	// fill histos hInvMassAll and hInvMassRes
	189	hInvMassAll->Fill(invMass);
	190	hInvMassRes->Fill(invMass);
8904b7e3	191	Float_t pxRes = Tracks_fPxRec[t1]+Tracks_fPxRec[t2];
	192	Float_t pyRes = Tracks_fPyRec[t1]+Tracks_fPyRec[t2];
	193	Float_t pzRes = Tracks_fPzRec[t1]+Tracks_fPzRec[t2];
	194	if (invMass > massMin && invMass < massMax) {
	195	EventInMass++;
	196	Float_t rapRes = 0;
	197	if (ResType == 553) {
	198	rapRes = rapParticle(pxRes,pyRes,pzRes,UpsilonMass);
	199	} else {
	200	rapRes = rapParticle(pxRes,pyRes,pzRes,JPsiMass);
	201	}
	202	hRapResonance->Fill(rapRes);
	203	hPtResonance->Fill(TMath::Sqrt(pxRespxRes+pyRespyRes));
	204	}
	205
c2a319d4	206	} //if ((Tracks_fCharge[t1] * Tracks_fCharge[t2]) == -1)
8904b7e3	207	} // if ((Tracks_fChi2[t2] < Chi2Cut) && (pt2 > PtCutMin) && (pt2 < PtCutMax))
c2a319d4	208	} // for (Int_t t2 = t1 + 1; t2 < Tracks_; t2++)
8904b7e3	209	} // if ((Tracks_fChi2[t1] < Chi2Cut) && (pt1 > PtCutMin)&& (pt1 < PtCutMax) )
c2a319d4	210	} // for (Int_t t1 = 0; t1 < Tracks_; t1++)
8904b7e3	211	hNumberOfTrack->Fill(NumberOfTrack);
c2a319d4	212	} // for (Int_t event = FirstEvent;
	213
	214	histoFile->Write();
	215	histoFile->Close();
8904b7e3	216
	217	cout << "MUONmassPlot " << endl;
	218	cout << "FirstEvent " << FirstEvent << endl;
	219	cout << "LastEvent " << LastEvent << endl;
	220	cout << "ResType " << ResType << endl;
	221	// cout << "Nsig " << Nsig << endl;
	222	cout << "Chi2Cut " << Chi2Cut << endl;
	223	cout << "PtCutMin " << PtCutMin << endl;
	224	cout << "PtCutMax " << PtCutMax << endl;
	225	cout << "massMin " << massMin << endl;
	226	cout << "massMax " << massMax << endl;
	227	cout << "EventInMass " << EventInMass << endl;
c2a319d4	228	}
c2a319d4	229
8904b7e3	230
8904b7e3	231	Float_t MuPlusMuMinusMass(Float_t Px1, Float_t Py1, Float_t Pz1, Float_t Px2, Float_t Py2, Float_t Pz2, Float_t muonMass)
c2a319d4	232	{
8904b7e3	233	// return invariant mass for particle 1 & 2 whose masses are equal to muonMass
8904b7e3	234
c2a319d4	235	Float_t e1 = TMath::Sqrt(muonMass * muonMass + Px1 * Px1 + Py1 * Py1 + Pz1 * Pz1);
	236	Float_t e2 = TMath::Sqrt(muonMass * muonMass + Px2 * Px2 + Py2 * Py2 + Pz2 * Pz2);
	237	return (TMath::Sqrt(2.0 * (muonMass * muonMass + e1 * e2 - Px1 * Px2 - Py1 * Py2 - Pz1 * Pz2)));
	238	}
8904b7e3	239
	240	Float_t rapParticle(Float_t Px, Float_t Py, Float_t Pz, Float_t Mass)
	241	{
	242	// return rapidity for particle
	243
	244	// Particle energy
	245	Float_t Energy = TMath::Sqrt(PxPx+PyPy+PzPz+MassMass);
	246	// Rapidity
	247	Float_t Rapidity = 0.5*TMath::Log((Energy+Pz)/(Energy-Pz));
	248	return Rapidity;
	249	}