[u/mrichter/AliRoot.git] / ANALYSIS / AliMuonAnalysis.cxx

#include "AliMuonAnalysis.h"
//________________________________
///////////////////////////////////////////////////////////
//
// class AliMuonAnalysis
//
// MUON Analysis
//
//
// 
// finck@subatech.in2p3.fr
//
///////////////////////////////////////////////////////////
/*********************************************************/

#include <TString.h>
#include <TParticle.h>

#include <AliStack.h>
#include <AliAOD.h>
#include <AliAODParticle.h>
#include <AliAODParticleCut.h>

#include <AliESDtrack.h>
#include <AliESD.h>

#include "TFile.h"
#include "TH1.h"
#include "TH2.h"

ClassImp(AliMuonAnalysis)

AliMuonAnalysis::AliMuonAnalysis():
  fHistoFile(0x0),
  fHPtMuon(0x0),
  fHPtMuonPlus(0x0),
  fHPtMuonMinus(0x0),
  fHPMuon(0x0),
  fHInvMassAll(0x0),
  fHRapMuon(0x0),
  fHRapResonance(0x0),
  fHPtResonance(0x0),
  fHInvMassAllvsPt(0x0),
  fPartCut(0x0)
{
 //ctor
}

/*********************************************************/
AliMuonAnalysis::~AliMuonAnalysis()
{
 //dtor
  delete fPartCut;
  delete fHistoFile;
  delete fHPtMuon;
  delete fHPtMuonPlus;
  delete fHPtMuonMinus;
  delete fHPMuon;
  delete fHInvMassAll;
  delete fHRapMuon;
  delete fHRapResonance;
  delete fHPtResonance;
  delete fHInvMassAllvsPt;
}
/*********************************************************/

Int_t AliMuonAnalysis::Init()
{
  //Initilizes analysis
  Info("Init","Histo initialized for MUON Analysis");

  fHistoFile         = new TFile("MUONmassPlot.root", "RECREATE");
  fHPtMuon           = new TH1F("hPtMuon", "Muon Pt (GeV/c)", 100, 0., 20.);
  fHPMuon            = new TH1F("hPMuon", "Muon P (GeV/c)", 100, 0., 200.);
  fHPtMuonPlus       = new TH1F("hPtMuonPlus", "Muon+ Pt (GeV/c)", 100, 0., 20.);
  fHPtMuonMinus      = new TH1F("hPtMuonMinus", "Muon- Pt (GeV/c)", 100, 0., 20.);
  fHInvMassAll       = new TH1F("hInvMassAll", "Mu+Mu- invariant mass (GeV/c2)", 480, 0., 12.);
  fHRapMuon          = new TH1F("hRapMuon"," Muon Rapidity",50,-4.5,-2);
  fHRapResonance     = new TH1F("hRapResonance"," Resonance Rapidity",50,-4.5,-2);
  fHPtResonance      = new TH1F("hPtResonance", "Resonance Pt (GeV/c)", 100, 0., 20.);
  fHInvMassAllvsPt = new TH2F("hInvMassAll_vs_Pt","hInvMassAll_vs_Pt",480,0.,12.,80,0.,20.);

  return 0; 
}
/*********************************************************/

Int_t AliMuonAnalysis::ProcessEvent(AliAOD* aodrec, AliAOD* aodsim)
{
  //
  // process the event
  // 
  if (aodrec) {
    GetInvMass(aodrec);
    //    Info("ProcessEvent","Inv Mass Rec");
  }  

  if (aodsim) {
    //     Info("ProcessEvent","aodsim not implemented");
  }  
  
  return 0;
  
}

/*********************************************************/

Int_t AliMuonAnalysis::Finish()
{
  //Finish analysis and writes results
  Info("Finish","Histo writing for MUON Analysis");

  fHistoFile->Write();
  fHistoFile->Close();

  return 0;
}
/*********************************************************/

void AliMuonAnalysis::GetInvMass(AliAOD* aod)
{
  // get the invariant mass distribution
  // from the oad events
 
  TLorentzVector lorV1, lorV2, lorVtot;
  Float_t massMin = 9.17;
  Float_t massMax = 9.77;
  Int_t charge1, charge2;

//returns flow parameters: v2 and event plane
  if (aod == 0x0) {
     Error("AliMuonAnalysis::GetInvMass","Pointer to AOD is NULL");
     return;
  }
   
  Int_t nPart = aod->GetNumberOfParticles();
  
  for (Int_t iPart1 = 0; iPart1 < nPart; iPart1++)  {
    AliAODParticle* aodPart1 =  (AliAODParticle*)aod->GetParticle(iPart1);

    if (aodPart1 == 0x0) {
      Error("AliMuonAnalysis::GetInvMass","Cannot get particle %d", iPart1);
      continue;
    }

    lorV1.SetPxPyPzE(aodPart1->Px(),
			 aodPart1->Py(),
			 aodPart1->Pz(),
			 aodPart1->E());

    fHPtMuon->Fill(lorV1.Pt());
    fHPMuon->Fill(lorV1.P());

    charge1 = TMath::Sign(1,aodPart1->GetPdgCode());

    if (charge1 > 0) {
      fHPtMuonPlus->Fill(lorV1.Pt());
    } else {
      fHPtMuonMinus->Fill(lorV1.Pt());
    }
    fHRapMuon->Fill(lorV1.Rapidity());
    for (Int_t iPart2 = iPart1 + 1; iPart2 < nPart; iPart2++)  {

      AliAODParticle* aodPart2 =  (AliAODParticle*)aod->GetParticle(iPart2);

      lorV2.SetPxPyPzE(aodPart2->Px(),
			 aodPart2->Py(),
			 aodPart2->Pz(),
			 aodPart2->E());

      charge2 = TMath::Sign(1,aodPart2->GetPdgCode());

      if ((charge1 * charge2) == -1) {

	lorVtot = lorV1;
	lorVtot += lorV2;
	Float_t invMass = lorVtot.M();

	fHInvMassAll->Fill(invMass);
	fHInvMassAllvsPt->Fill(invMass,lorVtot.Pt());

	if (invMass > massMin && invMass < massMax) {
	  fHRapResonance->Fill(lorVtot.Rapidity());
	  fHPtResonance->Fill(lorVtot.Pt());
	}
      }

    }
  }
}
Commit	Line	Data
beb1c41d	1	#include "AliMuonAnalysis.h"
	2	//________________________________
	3	///////////////////////////////////////////////////////////
	4	//
	5	// class AliMuonAnalysis
	6	//
	7	// MUON Analysis
	8	//
	9	//
	10	//
	11	// finck@subatech.in2p3.fr
	12	//
	13	///////////////////////////////////////////////////////////
	14	/*********************************************************/
	15
	16	#include <TString.h>
	17	#include <TParticle.h>
	18
	19	#include <AliStack.h>
	20	#include <AliAOD.h>
	21	#include <AliAODParticle.h>
	22	#include <AliAODParticleCut.h>
	23
	24	#include <AliESDtrack.h>
	25	#include <AliESD.h>
	26
	27	#include "TFile.h"
	28	#include "TH1.h"
	29	#include "TH2.h"
	30
	31	ClassImp(AliMuonAnalysis)
	32
	33	AliMuonAnalysis::AliMuonAnalysis():
62411bd0	34	fHistoFile(0x0),
	35	fHPtMuon(0x0),
	36	fHPtMuonPlus(0x0),
	37	fHPtMuonMinus(0x0),
	38	fHPMuon(0x0),
	39	fHInvMassAll(0x0),
	40	fHRapMuon(0x0),
	41	fHRapResonance(0x0),
	42	fHPtResonance(0x0),
	43	fHInvMassAllvsPt(0x0),
	44	fPartCut(0x0)
beb1c41d	45	{
	46	//ctor
	47	}
	48
	49	/*********************************************************/
	50	AliMuonAnalysis::~AliMuonAnalysis()
	51	{
	52	//dtor
	53	delete fPartCut;
	54	delete fHistoFile;
	55	delete fHPtMuon;
	56	delete fHPtMuonPlus;
	57	delete fHPtMuonMinus;
	58	delete fHPMuon;
	59	delete fHInvMassAll;
	60	delete fHRapMuon;
	61	delete fHRapResonance;
	62	delete fHPtResonance;
7255f192	63	delete fHInvMassAllvsPt;
beb1c41d	64	}
	65	/*********************************************************/
	66
	67	Int_t AliMuonAnalysis::Init()
	68	{
7255f192	69	//Initilizes analysis
beb1c41d	70	Info("Init","Histo initialized for MUON Analysis");
	71
	72	fHistoFile = new TFile("MUONmassPlot.root", "RECREATE");
	73	fHPtMuon = new TH1F("hPtMuon", "Muon Pt (GeV/c)", 100, 0., 20.);
	74	fHPMuon = new TH1F("hPMuon", "Muon P (GeV/c)", 100, 0., 200.);
	75	fHPtMuonPlus = new TH1F("hPtMuonPlus", "Muon+ Pt (GeV/c)", 100, 0., 20.);
	76	fHPtMuonMinus = new TH1F("hPtMuonMinus", "Muon- Pt (GeV/c)", 100, 0., 20.);
	77	fHInvMassAll = new TH1F("hInvMassAll", "Mu+Mu- invariant mass (GeV/c2)", 480, 0., 12.);
	78	fHRapMuon = new TH1F("hRapMuon"," Muon Rapidity",50,-4.5,-2);
	79	fHRapResonance = new TH1F("hRapResonance"," Resonance Rapidity",50,-4.5,-2);
	80	fHPtResonance = new TH1F("hPtResonance", "Resonance Pt (GeV/c)", 100, 0., 20.);
7255f192	81	fHInvMassAllvsPt = new TH2F("hInvMassAll_vs_Pt","hInvMassAll_vs_Pt",480,0.,12.,80,0.,20.);
beb1c41d	82
	83	return 0;
	84	}
	85	/*********************************************************/
	86
	87	Int_t AliMuonAnalysis::ProcessEvent(AliAOD* aodrec, AliAOD* aodsim)
	88	{
7255f192	89	//
	90	// process the event
	91	//
beb1c41d	92	if (aodrec) {
	93	GetInvMass(aodrec);
	94	// Info("ProcessEvent","Inv Mass Rec");
	95	}
	96
	97	if (aodsim) {
	98	// Info("ProcessEvent","aodsim not implemented");
	99	}
	100
	101	return 0;
	102
	103	}
	104
	105	/*********************************************************/
	106
	107	Int_t AliMuonAnalysis::Finish()
	108	{
	109	//Finish analysis and writes results
	110	Info("Finish","Histo writing for MUON Analysis");
	111
	112	fHistoFile->Write();
	113	fHistoFile->Close();
	114
	115	return 0;
	116	}
	117	/*********************************************************/
	118
	119	void AliMuonAnalysis::GetInvMass(AliAOD* aod)
	120	{
7255f192	121	// get the invariant mass distribution
	122	// from the oad events
	123
beb1c41d	124	TLorentzVector lorV1, lorV2, lorVtot;
	125	Float_t massMin = 9.17;
	126	Float_t massMax = 9.77;
	127	Int_t charge1, charge2;
	128
	129	//returns flow parameters: v2 and event plane
	130	if (aod == 0x0) {
	131	Error("AliMuonAnalysis::GetInvMass","Pointer to AOD is NULL");
	132	return;
	133	}
	134
	135	Int_t nPart = aod->GetNumberOfParticles();
	136
	137	for (Int_t iPart1 = 0; iPart1 < nPart; iPart1++) {
	138	AliAODParticle* aodPart1 = (AliAODParticle*)aod->GetParticle(iPart1);
	139
	140	if (aodPart1 == 0x0) {
	141	Error("AliMuonAnalysis::GetInvMass","Cannot get particle %d", iPart1);
	142	continue;
	143	}
	144
065e0e69	145	lorV1.SetPxPyPzE(aodPart1->Px(),
	146	aodPart1->Py(),
	147	aodPart1->Pz(),
	148	aodPart1->E());
beb1c41d	149
	150	fHPtMuon->Fill(lorV1.Pt());
	151	fHPMuon->Fill(lorV1.P());
	152
	153	charge1 = TMath::Sign(1,aodPart1->GetPdgCode());
	154
	155	if (charge1 > 0) {
	156	fHPtMuonPlus->Fill(lorV1.Pt());
	157	} else {
	158	fHPtMuonMinus->Fill(lorV1.Pt());
	159	}
	160	fHRapMuon->Fill(lorV1.Rapidity());
	161	for (Int_t iPart2 = iPart1 + 1; iPart2 < nPart; iPart2++) {
	162
	163	AliAODParticle* aodPart2 = (AliAODParticle*)aod->GetParticle(iPart2);
	164
065e0e69	165	lorV2.SetPxPyPzE(aodPart2->Px(),
	166	aodPart2->Py(),
	167	aodPart2->Pz(),
	168	aodPart2->E());
	169
beb1c41d	170	charge2 = TMath::Sign(1,aodPart2->GetPdgCode());
	171
	172	if ((charge1 * charge2) == -1) {
	173
065e0e69	174	lorVtot = lorV1;
065e0e69	175	lorVtot += lorV2;
beb1c41d	176	Float_t invMass = lorVtot.M();
	177
	178	fHInvMassAll->Fill(invMass);
7255f192	179	fHInvMassAllvsPt->Fill(invMass,lorVtot.Pt());
beb1c41d	180
	181	if (invMass > massMin && invMass < massMax) {
	182	fHRapResonance->Fill(lorVtot.Rapidity());
	183	fHPtResonance->Fill(lorVtot.Pt());
	184	}
	185	}
	186
	187	}
	188	}
	189	}