[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
//   NSigma (default 3)
//      the number of combinations is counted around the resonance mass
//      within +/- NSigma times the nominal sigma's
//      (0.099 GeV for Upsilon, 0.0615 GeV for J/Psi)
//   Chi2Cut (default 100)
//      to keep only tracks with chi2 per d.o.f. < Chi2Cut
//   PtCut (default 1)
//      to keep only tracks with transverse momentum > PtCut

// IMPORTANT NOTICE FOR USERS:
// under "root" or "root.exe", execute the following commands:
// 1. "gSystem->SetIncludePath("-I$ALICE_ROOT/MUON -I$ALICE_ROOT/STEER")" to get the right path at compilation time
// 2. ".x loadlibs.C" to load the shared libraries
// 3. ".L MUONrecoNtuple.C+"
// 4. ".x MUONmassPlot.C()" with the right arguments according to the list above

void MUONmassPlot(Int_t FirstEvent = 0, Int_t LastEvent = 0, Int_t ResType = 553, Float_t Nsig = 3., Float_t Chi2Cut = 100., Float_t PtCut = 1.)
{
  cout << "MUONmassPlot" << endl;
  cout << "FirstEvent" << FirstEvent << endl;
  cout << "LastEvent" << LastEvent << endl;
  cout << "ResType" << ResType << endl;
  cout << "Nsig" << Nsig << endl;
  cout << "Chi2Cut" << Chi2Cut << endl;
  cout << "PtCut" << PtCut << 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();
  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[5];
  Float_t         Tracks_fPxRec[5];
  Float_t         Tracks_fPyRec[5];
  Float_t         Tracks_fPzRec[5];
  Float_t         Tracks_fZRec[5];
  Float_t         Tracks_fZRec1[5];
  Int_t           Tracks_fNHits[5];
  Float_t         Tracks_fChi2[5];
  Float_t         Tracks_fPxGen[5];
  Float_t         Tracks_fPyGen[5];
  Float_t         Tracks_fPzGen[5];
  UInt_t          Tracks_fUniqueID[5];
  UInt_t          Tracks_fBits[5];

  //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 *hChi2PerDof = new TH1F("hChi2PerDof", "Muon track chi2/d.o.f.", 100, 0., 20.);
  TH1F *hInvMassAll = new TH1F("hInvMassAll", "Mu+Mu- invariant mass (GeV/c2)", 240, 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, 1., 5.);

  // Loop over events
  for (Int_t event = FirstEvent; event <= TMath::Min(LastEvent, nentries - 1); event++) {
    // get current event
    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]);
      // chi2 per d.o.f.
      Float_t ch1 = Tracks_fChi2[t1] / (2.0 * Tracks_fNHits[t1] - 5);
      // condition for good track (Chi2Cut and PtCut)
      if ((ch1 < Chi2Cut) && (pt1 > PtCut)) {
	// fill histos hPtMuon and hChi2PerDof
	hPtMuon->Fill(pt1);
	hChi2PerDof->Fill(ch1);
	// 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 > PtCut)) {
	    // 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]);
	      // fill histos hInvMassAll and hInvMassRes
	      hInvMassAll->Fill(invMass);
	      hInvMassRes->Fill(invMass);
	    } //if ((Tracks_fCharge[t1] * Tracks_fCharge[t2]) == -1)
	  } // if ((Tracks_fChi2[t2] < Chi2Cut) && (pt2 > PtCut))
	} // for (Int_t t2 = t1 + 1; t2 < Tracks_; t2++)
      } // if ((Tracks_fChi2[t1] < Chi2Cut) && (pt1 > PtCut))
    } // for (Int_t t1 = 0; t1 < Tracks_; t1++)
  } // for (Int_t event = FirstEvent;

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

Float_t MuPlusMuMinusMass(Float_t Px1, Float_t Py1, Float_t Pz1, Float_t Px2, Float_t Py2, Float_t Pz2)
{
  Float_t muonMass = 0.10566;
  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)));
}
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
	17	// NSigma (default 3)
	18	// the number of combinations is counted around the resonance mass
	19	// within +/- NSigma times the nominal sigma's
	20	// (0.099 GeV for Upsilon, 0.0615 GeV for J/Psi)
	21	// Chi2Cut (default 100)
	22	// to keep only tracks with chi2 per d.o.f. < Chi2Cut
	23	// PtCut (default 1)
	24	// to keep only tracks with transverse momentum > PtCut
	25
eff57104	26	// IMPORTANT NOTICE FOR USERS:
	27	// under "root" or "root.exe", execute the following commands:
	28	// 1. "gSystem->SetIncludePath("-I$ALICE_ROOT/MUON -I$ALICE_ROOT/STEER")" to get the right path at compilation time
	29	// 2. ".x loadlibs.C" to load the shared libraries
	30	// 3. ".L MUONrecoNtuple.C+"
	31	// 4. ".x MUONmassPlot.C()" with the right arguments according to the list above
	32
c2a319d4	33	void MUONmassPlot(Int_t FirstEvent = 0, Int_t LastEvent = 0, Int_t ResType = 553, Float_t Nsig = 3., Float_t Chi2Cut = 100., Float_t PtCut = 1.)
	34	{
	35	cout << "MUONmassPlot" << endl;
	36	cout << "FirstEvent" << FirstEvent << endl;
	37	cout << "LastEvent" << LastEvent << endl;
	38	cout << "ResType" << ResType << endl;
	39	cout << "Nsig" << Nsig << endl;
	40	cout << "Chi2Cut" << Chi2Cut << endl;
	41	cout << "PtCut" << PtCut << endl;
	42
	43	//////////////////////////////////////////////////////////
	44	// This file has been automatically generated
	45	// (Thu Sep 21 14:53:11 2000 by ROOT version2.25/02)
	46	// from TTree MUONtrackReco/MUONtrackReco
	47	// found on file: MUONtrackReco.root
	48	//////////////////////////////////////////////////////////
	49
	50
	51	//Reset ROOT and connect tree file
	52	gROOT->Reset();
	53	TFile f = (TFile)gROOT->GetListOfFiles()->FindObject("MUONtrackReco.root");
	54	if (!f) {
	55	f = new TFile("MUONtrackReco.root");
	56	}
	57	TTree MUONtrackReco = (TTree)gDirectory->Get("MUONtrackReco");
eff57104	58	MUONtrackReco->SetMakeClass(1);
c2a319d4	59
	60	//Declaration of leaves types
	61	Int_t fEvent;
	62	UInt_t fUniqueID;
	63	UInt_t fBits;
	64	Int_t Tracks_;
	65	Int_t Tracks_fCharge[5];
	66	Float_t Tracks_fPxRec[5];
	67	Float_t Tracks_fPyRec[5];
	68	Float_t Tracks_fPzRec[5];
	69	Float_t Tracks_fZRec[5];
	70	Float_t Tracks_fZRec1[5];
	71	Int_t Tracks_fNHits[5];
	72	Float_t Tracks_fChi2[5];
	73	Float_t Tracks_fPxGen[5];
	74	Float_t Tracks_fPyGen[5];
	75	Float_t Tracks_fPzGen[5];
	76	UInt_t Tracks_fUniqueID[5];
	77	UInt_t Tracks_fBits[5];
	78
	79	//Set branch addresses
	80	//MUONtrackReco->SetBranchAddress("Header",&Header);
	81	MUONtrackReco->SetBranchAddress("fEvent",&fEvent);
	82	MUONtrackReco->SetBranchAddress("fUniqueID",&fUniqueID);
	83	MUONtrackReco->SetBranchAddress("fBits",&fBits);
eff57104	84	// MUONtrackReco->SetBranchAddress("Tracks_",&Tracks_);
eff57104	85	MUONtrackReco->SetBranchAddress("Tracks",&Tracks_);
c2a319d4	86	MUONtrackReco->SetBranchAddress("Tracks.fCharge",Tracks_fCharge);
	87	MUONtrackReco->SetBranchAddress("Tracks.fPxRec",Tracks_fPxRec);
	88	MUONtrackReco->SetBranchAddress("Tracks.fPyRec",Tracks_fPyRec);
	89	MUONtrackReco->SetBranchAddress("Tracks.fPzRec",Tracks_fPzRec);
	90	MUONtrackReco->SetBranchAddress("Tracks.fZRec",Tracks_fZRec);
	91	MUONtrackReco->SetBranchAddress("Tracks.fZRec1",Tracks_fZRec1);
	92	MUONtrackReco->SetBranchAddress("Tracks.fNHits",Tracks_fNHits);
	93	MUONtrackReco->SetBranchAddress("Tracks.fChi2",Tracks_fChi2);
	94	MUONtrackReco->SetBranchAddress("Tracks.fPxGen",Tracks_fPxGen);
	95	MUONtrackReco->SetBranchAddress("Tracks.fPyGen",Tracks_fPyGen);
	96	MUONtrackReco->SetBranchAddress("Tracks.fPzGen",Tracks_fPzGen);
	97	MUONtrackReco->SetBranchAddress("Tracks.fUniqueID",Tracks_fUniqueID);
	98	MUONtrackReco->SetBranchAddress("Tracks.fBits",Tracks_fBits);
	99
	100	// This is the loop skeleton
	101	// To read only selected branches, Insert statements like:
	102	// MUONtrackReco->SetBranchStatus("*",0); // disable all branches
	103	// TTreePlayer->SetBranchStatus("branchname",1); // activate branchname
	104
	105	Int_t nentries = MUONtrackReco->GetEntries();
	106
	107	Int_t nbytes = 0;
	108	// for (Int_t i=0; i<nentries;i++) {
	109	// nbytes += MUONtrackReco->GetEntry(i);
	110	// }
	111
	112	/////////////////////////////////////////////////////////////////
	113	// Here comes the specialized part for MUONmassPlot
	114	/////////////////////////////////////////////////////////////////
	115
	116	// File for histograms and histogram booking
	117	TFile *histoFile = new TFile("MUONmassPlot.root", "RECREATE");
	118	TH1F *hPtMuon = new TH1F("hPtMuon", "Muon Pt (GeV/c)", 100, 0., 20.);
	119	TH1F *hChi2PerDof = new TH1F("hChi2PerDof", "Muon track chi2/d.o.f.", 100, 0., 20.);
	120	TH1F *hInvMassAll = new TH1F("hInvMassAll", "Mu+Mu- invariant mass (GeV/c2)", 240, 0., 12.);
	121	if (ResType = 553) TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around Upsilon", 60, 8., 11.);
	122	else TH1F *hInvMassRes = new TH1F("hInvMassRes", "Mu+Mu- invariant mass (GeV/c2) around J/Psi", 80, 1., 5.);
	123
	124	// Loop over events
	125	for (Int_t event = FirstEvent; event <= TMath::Min(LastEvent, nentries - 1); event++) {
	126	// get current event
	127	nbytes += MUONtrackReco->GetEntry(event);
	128	// loop over all reconstructed tracks (also first track of combination)
	129	for (Int_t t1 = 0; t1 < Tracks_; t1++) {
	130	// transverse momentum
	131	Float_t pt1 = TMath::Sqrt(Tracks_fPxRec[t1] * Tracks_fPxRec[t1] + Tracks_fPyRec[t1] * Tracks_fPyRec[t1]);
	132	// chi2 per d.o.f.
	133	Float_t ch1 = Tracks_fChi2[t1] / (2.0 * Tracks_fNHits[t1] - 5);
	134	// condition for good track (Chi2Cut and PtCut)
	135	if ((ch1 < Chi2Cut) && (pt1 > PtCut)) {
	136	// fill histos hPtMuon and hChi2PerDof
	137	hPtMuon->Fill(pt1);
	138	hChi2PerDof->Fill(ch1);
	139	// loop over second track of combination
	140	for (Int_t t2 = t1 + 1; t2 < Tracks_; t2++) {
	141	// transverse momentum
	142	Float_t pt2 = TMath::Sqrt(Tracks_fPxRec[t2] * Tracks_fPxRec[t2] + Tracks_fPyRec[t2] * Tracks_fPyRec[t2]);
	143	// chi2 per d.o.f.
	144	Float_t ch2 = Tracks_fChi2[t2] / (2.0 * Tracks_fNHits[t2] - 5);
	145	// condition for good track (Chi2Cut and PtCut)
	146	if ((ch2 < Chi2Cut) && (pt2 > PtCut)) {
	147	// condition for opposite charges
	148	if ((Tracks_fCharge[t1] * Tracks_fCharge[t2]) == -1) {
	149	// invariant mass
150	Float_t invMass = MuPlusMuMinusMass(Tracks_fPxRec[t1], Tracks_fPyRec[t1], Tracks_fPzRec[t1], Tracks_fPxRec[t2], Tracks_fPyRec[t2], Tracks_fPzRec[t2]);
151	// fill histos hInvMassAll and hInvMassRes
152	hInvMassAll->Fill(invMass);
153	hInvMassRes->Fill(invMass);
154	} //if ((Tracks_fCharge[t1] * Tracks_fCharge[t2]) == -1)
155	} // if ((Tracks_fChi2[t2] < Chi2Cut) && (pt2 > PtCut))
156	} // for (Int_t t2 = t1 + 1; t2 < Tracks_; t2++)
157	} // if ((Tracks_fChi2[t1] < Chi2Cut) && (pt1 > PtCut))
158	} // for (Int_t t1 = 0; t1 < Tracks_; t1++)
159	} // for (Int_t event = FirstEvent;
160
161	histoFile->Write();
162	histoFile->Close();
163	}
164
165	Float_t MuPlusMuMinusMass(Float_t Px1, Float_t Py1, Float_t Pz1, Float_t Px2, Float_t Py2, Float_t Pz2)
166	{
167	Float_t muonMass = 0.10566;
168	Float_t e1 = TMath::Sqrt(muonMass * muonMass + Px1 * Px1 + Py1 * Py1 + Pz1 * Pz1);
169	Float_t e2 = TMath::Sqrt(muonMass * muonMass + Px2 * Px2 + Py2 * Py2 + Pz2 * Pz2);
170	return (TMath::Sqrt(2.0 * (muonMass * muonMass + e1 * e2 - Px1 * Px2 - Py1 * Py2 - Pz1 * Pz2)));
171	}