[u/mrichter/AliRoot.git] / MFT / AliAnalysisTaskMFTExample.cxx

#include "TH1D.h"
#include "AliAODEvent.h"
#include "AliAODMCParticle.h"
#include "AliAODMCHeader.h"
#include "AliAnalysisTaskMFTExample.h"
#include "TDatabasePDG.h"
#include "AliMUONTrackParam.h"
#include "AliMUONTrackExtrap.h"
#include "AliAODDimuon.h"
#include "AliAODTrack.h"
#include "AliAODHeader.h"
#include "TClonesArray.h"
#include "AliMFTConstants.h"
#include "AliMFTSupport.h"
#include "TRandom.h"

ClassImp(AliAnalysisTaskMFTExample)

//====================================================================================================================================================

AliAnalysisTaskMFTExample::AliAnalysisTaskMFTExample() : 
  AliAnalysisTaskSE(),
  fVertexMode(0),
  fHistPtSingleMuons(0),
  fHistPtSingleMuonsFromJpsi(0),
  fHistPtDimuonsOS(0),
  fHistMassDimuonsOS(0),
  fHistPtDimuonsJpsi(0),
  fHistMassDimuonsJpsi(0),
  fHistResidualXVtxJpsi(0),
  fHistResidualYVtxJpsi(0),
  fHistResidualZVtxJpsi(0) {
  
  // Default constructor

  for (Int_t i=0; i<3; i++) fPrimaryVertex[i] = 0;
  fVtxResolutionITS[0] = 5.e-4;
  fVtxResolutionITS[1] = 5.e-4;
  fVtxResolutionITS[2] = 4.e-4;

}

//====================================================================================================================================================

AliAnalysisTaskMFTExample::AliAnalysisTaskMFTExample(const Char_t *name) : 
  AliAnalysisTaskSE(name),
  fVertexMode(0),
  fHistPtSingleMuons(0),
  fHistPtSingleMuonsFromJpsi(0),
  fHistPtDimuonsOS(0),
  fHistMassDimuonsOS(0),
  fHistPtDimuonsJpsi(0),
  fHistMassDimuonsJpsi(0),
  fHistResidualXVtxJpsi(0),
  fHistResidualYVtxJpsi(0),
  fHistResidualZVtxJpsi(0) {

  // Constructor

  for (Int_t i=0; i<3; i++) fPrimaryVertex[i] = 0;
  fVtxResolutionITS[0] = 5.e-4;
  fVtxResolutionITS[1] = 5.e-4;
  fVtxResolutionITS[2] = 4.e-4;

  // Define input and output slots here
  DefineOutput(1, TH1D::Class());
  DefineOutput(2, TH1D::Class());
  DefineOutput(3, TH1D::Class());
  DefineOutput(4, TH1D::Class());
  DefineOutput(5, TH1D::Class());
  DefineOutput(6, TH1D::Class());
  DefineOutput(7, TH1D::Class());
  DefineOutput(8, TH1D::Class());
  DefineOutput(9, TH1D::Class());

}

//====================================================================================================================================================

void AliAnalysisTaskMFTExample::UserCreateOutputObjects() {

  // Called once

  fHistPtSingleMuons = new TH1D("fHistPtSingleMuons","p_{T} of single muons (All)", 100, 0, 10);
  fHistPtSingleMuons -> SetXTitle("p_{T}  [GeV/c]");
  fHistPtSingleMuons -> Sumw2();

  fHistPtSingleMuonsFromJpsi = new TH1D("fHistPtSingleMuonsFromJpsi", "p_{T} of single muons (from J/#psi)", 100, 0, 10);
  fHistPtSingleMuonsFromJpsi -> SetXTitle("p_{T}  [GeV/c]");
  fHistPtSingleMuonsFromJpsi -> Sumw2();

  fHistMassDimuonsOS = new TH1D("fHistMassDimuonsOS", "Mass of OS dimuons", 500, 0, 10);
  fHistMassDimuonsOS -> SetXTitle("Mass  [GeV/c^{2}]");
  fHistMassDimuonsOS -> Sumw2();

  fHistMassDimuonsJpsi = new TH1D("fHistMassDimuonsJpsi", "Mass of J/#psi dimuons", 500, 0, 10);
  fHistMassDimuonsJpsi -> SetXTitle("Mass  [GeV/c^{2}]");
  fHistMassDimuonsJpsi -> Sumw2();

  fHistPtDimuonsOS = new TH1D("fHistPtDimuonsOS", "p_{T} of OS dimuons", 100, 0, 10);
  fHistPtDimuonsOS -> SetXTitle("p_{T}  [GeV/c]");
  fHistPtDimuonsOS -> Sumw2();

  fHistPtDimuonsJpsi = new TH1D("fHistPtDimuonsJpsi", "p_{T} of J/#psi dimuons", 100, 0, 10);
  fHistPtDimuonsJpsi -> SetXTitle("p_{T}  [GeV/c]");
  fHistPtDimuonsJpsi -> Sumw2();

  fHistResidualXVtxJpsi = new TH1D("fHistResidualXVtxJpsi", "J/#psi vertex residuals along x", 100, -100, 100);
  fHistResidualXVtxJpsi -> SetXTitle("Residual  [#mum]");
  fHistResidualXVtxJpsi -> Sumw2();

  fHistResidualYVtxJpsi = new TH1D("fHistResidualYVtxJpsi", "J/#psi vertex residuals along y", 100, -100, 100);
  fHistResidualYVtxJpsi -> SetXTitle("Residual  [#mum]");
  fHistResidualYVtxJpsi -> Sumw2();

  fHistResidualZVtxJpsi = new TH1D("fHistResidualZVtxJpsi", "J/#psi vertex residuals along z", 100, -1000, 1000);
  fHistResidualZVtxJpsi -> SetXTitle("Residual  [#mum]");
  fHistResidualZVtxJpsi -> Sumw2();

  PostData(1, fHistPtSingleMuons);
  PostData(2, fHistPtSingleMuonsFromJpsi);
  PostData(3, fHistMassDimuonsOS);
  PostData(4, fHistMassDimuonsJpsi);
  PostData(5, fHistPtDimuonsOS);
  PostData(6, fHistPtDimuonsJpsi);
  PostData(7, fHistResidualXVtxJpsi);
  PostData(8, fHistResidualYVtxJpsi);
  PostData(9, fHistResidualZVtxJpsi);
  
}

//====================================================================================================================================================

void AliAnalysisTaskMFTExample::UserExec(Option_t *) {

  // Main loop
  // Called for each event	

  AliAODEvent *aodEv = dynamic_cast<AliAODEvent*>(InputEvent());  
  if (!aodEv) return;

  aodEv->GetHeader()->InitMagneticField();
  AliMUONTrackExtrap::SetField();

  TClonesArray *stackMC = (TClonesArray*) (aodEv->GetList()->FindObject(AliAODMCParticle::StdBranchName()));

  // Getting primary vertex, either from the generation or from the reconstruction -------------------

  if (fVertexMode == kGenerated) {
    AliAODMCHeader *mcHeader = (AliAODMCHeader*) (aodEv->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
    mcHeader->GetVertex(fPrimaryVertex);
    for (Int_t i=0; i<3; i++) fPrimaryVertex[i] = gRandom->Gaus(fPrimaryVertex[i], fVtxResolutionITS[i]);
  }
  else if (fVertexMode == kReconstructed) {
    aodEv->GetPrimaryVertex()->GetXYZ(fPrimaryVertex);
  }

  // -------------------------------------------------------------------------------------------------

  AliAODTrack *recMuon1=0, *recMuon2=0;
  AliAODMCParticle *mcMuon1=0, *mcMuon2=0;
  AliAODMCParticle *mcMother1=0;
  Bool_t isMuon1FromJpsi=0;

  // Loop over MUON+MFT muons

  //--------------------------------------------------------------------------------------------------

  for (Int_t iTrack=0; iTrack<aodEv->GetNumberOfTracks(); iTrack++) { 

    if (!(aodEv->GetTrack(iTrack)->IsMuonGlobalTrack())) continue;

    recMuon1 = aodEv->GetTrack(iTrack);
    isMuon1FromJpsi = kFALSE;
    
    // pt all muons
    fHistPtSingleMuons -> Fill(recMuon1->Pt());

    // pt muons from J/psi
    if (recMuon1->GetLabel() >= 0) {
      mcMuon1 = (AliAODMCParticle*) stackMC->At(recMuon1->GetLabel());
      if (mcMuon1) {
	if (mcMuon1->GetMother()>=0) {
	  mcMother1 = (AliAODMCParticle*) stackMC->At(mcMuon1->GetMother());
	  if (mcMother1->PdgCode()==443) {
	    isMuon1FromJpsi = kTRUE;
	    fHistPtSingleMuonsFromJpsi -> Fill(recMuon1->Pt());
	  }
	}
      }
    }

    for (Int_t jTrack=0; jTrack<iTrack; jTrack++) { 

      if (!(aodEv->GetTrack(jTrack)->IsMuonGlobalTrack())) continue;
      
      recMuon2 = aodEv->GetTrack(jTrack);
      
      AliAODDimuon *dimuon = new AliAODDimuon(recMuon1, recMuon2);
      if (dimuon->Charge()) continue;

      // pt and mass all OS dimuons
      fHistPtDimuonsOS   -> Fill(dimuon->Pt());
      fHistMassDimuonsOS -> Fill(dimuon->Mass());

      // pt and mass J/psi dimuons
      if (!isMuon1FromJpsi) continue;
      if (recMuon2->GetLabel() >= 0) {
	mcMuon2 = (AliAODMCParticle*) stackMC->At(recMuon2->GetLabel());
	if (mcMuon2) {
	  if (mcMuon2->GetMother() == mcMuon1->GetMother()) {
	    Double_t vertex[3]={0};
	    TLorentzVector kinem(0,0,0,0);
	    if (!AliMFTSupport::RefitAODDimuonWithCommonVertex(dimuon, vertex, kinem)) continue;
	    fHistPtDimuonsJpsi    -> Fill(kinem.Pt());
	    fHistMassDimuonsJpsi  -> Fill(kinem.M());
	    fHistResidualXVtxJpsi -> Fill(1.e4*(vertex[0] - fPrimaryVertex[0]));
	    fHistResidualYVtxJpsi -> Fill(1.e4*(vertex[1] - fPrimaryVertex[1]));
	    fHistResidualZVtxJpsi -> Fill(1.e4*(vertex[2] - fPrimaryVertex[2]));
	  }
	}
      }
   
    }   // end of loop on 2nd muon

  }   // end of loop on 1st muon
 
  //--------------------------------------------------------------------------------------------------
   
  PostData(1, fHistPtSingleMuons);
  PostData(2, fHistPtSingleMuonsFromJpsi);
  PostData(3, fHistMassDimuonsOS);
  PostData(4, fHistMassDimuonsJpsi);
  PostData(5, fHistPtDimuonsOS);
  PostData(6, fHistPtDimuonsJpsi);
  PostData(7, fHistResidualXVtxJpsi);
  PostData(8, fHistResidualYVtxJpsi);
  PostData(9, fHistResidualZVtxJpsi);

}

//====================================================================================================================================================

void AliAnalysisTaskMFTExample::Terminate(Option_t *) {

  // Draw result to the screen
  // Called once at the end of the query
  
}

//====================================================================================================================================================
Commit	Line	Data
d30a0b7b	1	#include "TH1D.h"
	2	#include "AliAODEvent.h"
	3	#include "AliAODMCParticle.h"
	4	#include "AliAODMCHeader.h"
	5	#include "AliAnalysisTaskMFTExample.h"
	6	#include "TDatabasePDG.h"
	7	#include "AliMUONTrackParam.h"
	8	#include "AliMUONTrackExtrap.h"
	9	#include "AliAODDimuon.h"
	10	#include "AliAODTrack.h"
	11	#include "AliAODHeader.h"
	12	#include "TClonesArray.h"
	13	#include "AliMFTConstants.h"
	14	#include "AliMFTSupport.h"
	15	#include "TRandom.h"
	16
	17	ClassImp(AliAnalysisTaskMFTExample)
	18
	19	//====================================================================================================================================================
	20
	21	AliAnalysisTaskMFTExample::AliAnalysisTaskMFTExample() :
	22	AliAnalysisTaskSE(),
	23	fVertexMode(0),
	24	fHistPtSingleMuons(0),
	25	fHistPtSingleMuonsFromJpsi(0),
	26	fHistPtDimuonsOS(0),
	27	fHistMassDimuonsOS(0),
	28	fHistPtDimuonsJpsi(0),
	29	fHistMassDimuonsJpsi(0),
	30	fHistResidualXVtxJpsi(0),
	31	fHistResidualYVtxJpsi(0),
	32	fHistResidualZVtxJpsi(0) {
	33
	34	// Default constructor
	35
	36	for (Int_t i=0; i<3; i++) fPrimaryVertex[i] = 0;
	37	fVtxResolutionITS[0] = 5.e-4;
	38	fVtxResolutionITS[1] = 5.e-4;
	39	fVtxResolutionITS[2] = 4.e-4;
	40
	41	}
	42
	43	//====================================================================================================================================================
	44
	45	AliAnalysisTaskMFTExample::AliAnalysisTaskMFTExample(const Char_t *name) :
	46	AliAnalysisTaskSE(name),
	47	fVertexMode(0),
	48	fHistPtSingleMuons(0),
	49	fHistPtSingleMuonsFromJpsi(0),
	50	fHistPtDimuonsOS(0),
	51	fHistMassDimuonsOS(0),
	52	fHistPtDimuonsJpsi(0),
	53	fHistMassDimuonsJpsi(0),
	54	fHistResidualXVtxJpsi(0),
	55	fHistResidualYVtxJpsi(0),
	56	fHistResidualZVtxJpsi(0) {
	57
	58	// Constructor
	59
	60	for (Int_t i=0; i<3; i++) fPrimaryVertex[i] = 0;
	61	fVtxResolutionITS[0] = 5.e-4;
	62	fVtxResolutionITS[1] = 5.e-4;
	63	fVtxResolutionITS[2] = 4.e-4;
	64
65	// Define input and output slots here
66	DefineOutput(1, TH1D::Class());
67	DefineOutput(2, TH1D::Class());
68	DefineOutput(3, TH1D::Class());
69	DefineOutput(4, TH1D::Class());
70	DefineOutput(5, TH1D::Class());
71	DefineOutput(6, TH1D::Class());
72	DefineOutput(7, TH1D::Class());
73	DefineOutput(8, TH1D::Class());
74	DefineOutput(9, TH1D::Class());
75
76	}
77
78	//====================================================================================================================================================
79
80	void AliAnalysisTaskMFTExample::UserCreateOutputObjects() {
81
82	// Called once
83
84	fHistPtSingleMuons = new TH1D("fHistPtSingleMuons","p_{T} of single muons (All)", 100, 0, 10);
85	fHistPtSingleMuons -> SetXTitle("p_{T} [GeV/c]");
86	fHistPtSingleMuons -> Sumw2();
87
88	fHistPtSingleMuonsFromJpsi = new TH1D("fHistPtSingleMuonsFromJpsi", "p_{T} of single muons (from J/#psi)", 100, 0, 10);
89	fHistPtSingleMuonsFromJpsi -> SetXTitle("p_{T} [GeV/c]");
90	fHistPtSingleMuonsFromJpsi -> Sumw2();
91
92	fHistMassDimuonsOS = new TH1D("fHistMassDimuonsOS", "Mass of OS dimuons", 500, 0, 10);
93	fHistMassDimuonsOS -> SetXTitle("Mass [GeV/c^{2}]");
94	fHistMassDimuonsOS -> Sumw2();
95
96	fHistMassDimuonsJpsi = new TH1D("fHistMassDimuonsJpsi", "Mass of J/#psi dimuons", 500, 0, 10);
97	fHistMassDimuonsJpsi -> SetXTitle("Mass [GeV/c^{2}]");
98	fHistMassDimuonsJpsi -> Sumw2();
99
100	fHistPtDimuonsOS = new TH1D("fHistPtDimuonsOS", "p_{T} of OS dimuons", 100, 0, 10);
101	fHistPtDimuonsOS -> SetXTitle("p_{T} [GeV/c]");
102	fHistPtDimuonsOS -> Sumw2();
103
104	fHistPtDimuonsJpsi = new TH1D("fHistPtDimuonsJpsi", "p_{T} of J/#psi dimuons", 100, 0, 10);
105	fHistPtDimuonsJpsi -> SetXTitle("p_{T} [GeV/c]");
106	fHistPtDimuonsJpsi -> Sumw2();
107
108	fHistResidualXVtxJpsi = new TH1D("fHistResidualXVtxJpsi", "J/#psi vertex residuals along x", 100, -100, 100);
109	fHistResidualXVtxJpsi -> SetXTitle("Residual [#mum]");
110	fHistResidualXVtxJpsi -> Sumw2();
111
112	fHistResidualYVtxJpsi = new TH1D("fHistResidualYVtxJpsi", "J/#psi vertex residuals along y", 100, -100, 100);
113	fHistResidualYVtxJpsi -> SetXTitle("Residual [#mum]");
114	fHistResidualYVtxJpsi -> Sumw2();
115
116	fHistResidualZVtxJpsi = new TH1D("fHistResidualZVtxJpsi", "J/#psi vertex residuals along z", 100, -1000, 1000);
117	fHistResidualZVtxJpsi -> SetXTitle("Residual [#mum]");
118	fHistResidualZVtxJpsi -> Sumw2();
119
120	PostData(1, fHistPtSingleMuons);
121	PostData(2, fHistPtSingleMuonsFromJpsi);
122	PostData(3, fHistMassDimuonsOS);
123	PostData(4, fHistMassDimuonsJpsi);
124	PostData(5, fHistPtDimuonsOS);
125	PostData(6, fHistPtDimuonsJpsi);
126	PostData(7, fHistResidualXVtxJpsi);
127	PostData(8, fHistResidualYVtxJpsi);
128	PostData(9, fHistResidualZVtxJpsi);
129
130	}
131
132	//====================================================================================================================================================
133
134	void AliAnalysisTaskMFTExample::UserExec(Option_t *) {
135
136	// Main loop
137	// Called for each event
138
139	AliAODEvent aodEv = dynamic_cast<AliAODEvent>(InputEvent());
140	if (!aodEv) return;
141
142	aodEv->GetHeader()->InitMagneticField();
143	AliMUONTrackExtrap::SetField();
144
145	TClonesArray stackMC = (TClonesArray) (aodEv->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
146
147	// Getting primary vertex, either from the generation or from the reconstruction -------------------
148
149	if (fVertexMode == kGenerated) {
150	AliAODMCHeader mcHeader = (AliAODMCHeader) (aodEv->GetList()->FindObject(AliAODMCHeader::StdBranchName()));
151	mcHeader->GetVertex(fPrimaryVertex);
152	for (Int_t i=0; i<3; i++) fPrimaryVertex[i] = gRandom->Gaus(fPrimaryVertex[i], fVtxResolutionITS[i]);
153	}
154	else if (fVertexMode == kReconstructed) {
155	aodEv->GetPrimaryVertex()->GetXYZ(fPrimaryVertex);
156	}
157
158	// -------------------------------------------------------------------------------------------------
159
160	AliAODTrack recMuon1=0, recMuon2=0;
161	AliAODMCParticle mcMuon1=0, mcMuon2=0;
162	AliAODMCParticle *mcMother1=0;
163	Bool_t isMuon1FromJpsi=0;
164
165	// Loop over MUON+MFT muons
166
167	//--------------------------------------------------------------------------------------------------
168
169	for (Int_t iTrack=0; iTrack<aodEv->GetNumberOfTracks(); iTrack++) {
170
171	if (!(aodEv->GetTrack(iTrack)->IsMuonGlobalTrack())) continue;
172
173	recMuon1 = aodEv->GetTrack(iTrack);
174	isMuon1FromJpsi = kFALSE;
175
176	// pt all muons
177	fHistPtSingleMuons -> Fill(recMuon1->Pt());
178
179	// pt muons from J/psi
180	if (recMuon1->GetLabel() >= 0) {
181	mcMuon1 = (AliAODMCParticle*) stackMC->At(recMuon1->GetLabel());
182	if (mcMuon1) {
183	if (mcMuon1->GetMother()>=0) {
184	mcMother1 = (AliAODMCParticle*) stackMC->At(mcMuon1->GetMother());
185	if (mcMother1->PdgCode()==443) {
186	isMuon1FromJpsi = kTRUE;
187	fHistPtSingleMuonsFromJpsi -> Fill(recMuon1->Pt());
188	}
189	}
190	}
191	}
192
193	for (Int_t jTrack=0; jTrack<iTrack; jTrack++) {
194
195	if (!(aodEv->GetTrack(jTrack)->IsMuonGlobalTrack())) continue;
196
197	recMuon2 = aodEv->GetTrack(jTrack);
198
199	AliAODDimuon *dimuon = new AliAODDimuon(recMuon1, recMuon2);
200	if (dimuon->Charge()) continue;
201
202	// pt and mass all OS dimuons
203	fHistPtDimuonsOS -> Fill(dimuon->Pt());
204	fHistMassDimuonsOS -> Fill(dimuon->Mass());
205
206	// pt and mass J/psi dimuons
207	if (!isMuon1FromJpsi) continue;
208	if (recMuon2->GetLabel() >= 0) {
209	mcMuon2 = (AliAODMCParticle*) stackMC->At(recMuon2->GetLabel());
210	if (mcMuon2) {
211	if (mcMuon2->GetMother() == mcMuon1->GetMother()) {
212	Double_t vertex[3]={0};
213	TLorentzVector kinem(0,0,0,0);
214	if (!AliMFTSupport::RefitAODDimuonWithCommonVertex(dimuon, vertex, kinem)) continue;
215	fHistPtDimuonsJpsi -> Fill(kinem.Pt());
216	fHistMassDimuonsJpsi -> Fill(kinem.M());
217	fHistResidualXVtxJpsi -> Fill(1.e4*(vertex[0] - fPrimaryVertex[0]));
218	fHistResidualYVtxJpsi -> Fill(1.e4*(vertex[1] - fPrimaryVertex[1]));
219	fHistResidualZVtxJpsi -> Fill(1.e4*(vertex[2] - fPrimaryVertex[2]));
220	}
221	}
222	}
223
224	} // end of loop on 2nd muon
225
226	} // end of loop on 1st muon
227
228	//--------------------------------------------------------------------------------------------------
229
230	PostData(1, fHistPtSingleMuons);
231	PostData(2, fHistPtSingleMuonsFromJpsi);
232	PostData(3, fHistMassDimuonsOS);
233	PostData(4, fHistMassDimuonsJpsi);
234	PostData(5, fHistPtDimuonsOS);
235	PostData(6, fHistPtDimuonsJpsi);
236	PostData(7, fHistResidualXVtxJpsi);
237	PostData(8, fHistResidualYVtxJpsi);
238	PostData(9, fHistResidualZVtxJpsi);
239
240	}
241
242	//====================================================================================================================================================
243
244	void AliAnalysisTaskMFTExample::Terminate(Option_t *) {
245
246	// Draw result to the screen
247	// Called once at the end of the query
248
249	}
250
251	//====================================================================================================================================================