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