[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 "AliMFTAnalysisTools.h"
#include "TRandom.h"
#include "TList.h"

ClassImp(AliAnalysisTaskMFTExample)

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

AliAnalysisTaskMFTExample::AliAnalysisTaskMFTExample() : 
  AliAnalysisTaskSE(),
  fVertexMode(0),
  fHistogramList(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),
  fHistogramList(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, TList::Class());

}

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

void AliAnalysisTaskMFTExample::UserCreateOutputObjects() {

  // Called once

  fHistogramList = new TList();
  fHistogramList->SetOwner(kTRUE);
  
  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();

  fHistogramList->Add(fHistPtSingleMuons);
  fHistogramList->Add(fHistPtSingleMuonsFromJpsi);
  fHistogramList->Add(fHistMassDimuonsOS);
  fHistogramList->Add(fHistMassDimuonsJpsi);
  fHistogramList->Add(fHistPtDimuonsOS);
  fHistogramList->Add(fHistPtDimuonsJpsi);
  fHistogramList->Add(fHistResidualXVtxJpsi);
  fHistogramList->Add(fHistResidualYVtxJpsi);
  fHistogramList->Add(fHistResidualZVtxJpsi);

  PostData(1, fHistogramList);

}

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

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()) {
	    AliAODDimuon *dimuon = new AliAODDimuon;
	    dimuon->SetMuons(recMuon1,recMuon2);
	    Double_t pca[3]={0};
	    Double_t pcaQuality=0;
	    TLorentzVector kinem(0,0,0,0);
	    if (!AliMFTAnalysisTools::CalculatePCA(dimuon, pca, pcaQuality, kinem)) continue;
	    fHistPtDimuonsJpsi    -> Fill(kinem.Pt());
	    fHistMassDimuonsJpsi  -> Fill(kinem.M());
	    fHistResidualXVtxJpsi -> Fill(1.e4*(pca[0] - fPrimaryVertex[0]));
	    fHistResidualYVtxJpsi -> Fill(1.e4*(pca[1] - fPrimaryVertex[1]));
	    fHistResidualZVtxJpsi -> Fill(1.e4*(pca[2] - fPrimaryVertex[2]));
	    delete dimuon;
	  }
	}
      }
   
    }   // end of loop on 2nd muon

  }   // end of loop on 1st muon
 
  //--------------------------------------------------------------------------------------------------
   
  PostData(1, fHistogramList);

}

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

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"
bffc7f8c	14	#include "AliMFTAnalysisTools.h"
d30a0b7b	15	#include "TRandom.h"
2cf4030c	16	#include "TList.h"
d30a0b7b	17
	18	ClassImp(AliAnalysisTaskMFTExample)
	19
	20	//====================================================================================================================================================
	21
	22	AliAnalysisTaskMFTExample::AliAnalysisTaskMFTExample() :
	23	AliAnalysisTaskSE(),
	24	fVertexMode(0),
2cf4030c	25	fHistogramList(0),
d30a0b7b	26	fHistPtSingleMuons(0),
	27	fHistPtSingleMuonsFromJpsi(0),
	28	fHistPtDimuonsOS(0),
	29	fHistMassDimuonsOS(0),
	30	fHistPtDimuonsJpsi(0),
	31	fHistMassDimuonsJpsi(0),
	32	fHistResidualXVtxJpsi(0),
	33	fHistResidualYVtxJpsi(0),
	34	fHistResidualZVtxJpsi(0) {
	35
	36	// Default constructor
	37
	38	for (Int_t i=0; i<3; i++) fPrimaryVertex[i] = 0;
	39	fVtxResolutionITS[0] = 5.e-4;
	40	fVtxResolutionITS[1] = 5.e-4;
	41	fVtxResolutionITS[2] = 4.e-4;
	42
	43	}
	44
	45	//====================================================================================================================================================
	46
	47	AliAnalysisTaskMFTExample::AliAnalysisTaskMFTExample(const Char_t *name) :
	48	AliAnalysisTaskSE(name),
	49	fVertexMode(0),
2cf4030c	50	fHistogramList(0),
d30a0b7b	51	fHistPtSingleMuons(0),
	52	fHistPtSingleMuonsFromJpsi(0),
	53	fHistPtDimuonsOS(0),
	54	fHistMassDimuonsOS(0),
	55	fHistPtDimuonsJpsi(0),
	56	fHistMassDimuonsJpsi(0),
	57	fHistResidualXVtxJpsi(0),
	58	fHistResidualYVtxJpsi(0),
	59	fHistResidualZVtxJpsi(0) {
	60
	61	// Constructor
	62
	63	for (Int_t i=0; i<3; i++) fPrimaryVertex[i] = 0;
	64	fVtxResolutionITS[0] = 5.e-4;
	65	fVtxResolutionITS[1] = 5.e-4;
	66	fVtxResolutionITS[2] = 4.e-4;
	67
	68	// Define input and output slots here
2cf4030c	69	DefineOutput(1, TList::Class());
d30a0b7b	70
	71	}
	72
	73	//====================================================================================================================================================
	74
	75	void AliAnalysisTaskMFTExample::UserCreateOutputObjects() {
	76
	77	// Called once
	78
2cf4030c	79	fHistogramList = new TList();
	80	fHistogramList->SetOwner(kTRUE);
	81
d30a0b7b	82	fHistPtSingleMuons = new TH1D("fHistPtSingleMuons","p_{T} of single muons (All)", 100, 0, 10);
	83	fHistPtSingleMuons -> SetXTitle("p_{T} [GeV/c]");
	84	fHistPtSingleMuons -> Sumw2();
	85
	86	fHistPtSingleMuonsFromJpsi = new TH1D("fHistPtSingleMuonsFromJpsi", "p_{T} of single muons (from J/#psi)", 100, 0, 10);
	87	fHistPtSingleMuonsFromJpsi -> SetXTitle("p_{T} [GeV/c]");
	88	fHistPtSingleMuonsFromJpsi -> Sumw2();
	89
	90	fHistMassDimuonsOS = new TH1D("fHistMassDimuonsOS", "Mass of OS dimuons", 500, 0, 10);
	91	fHistMassDimuonsOS -> SetXTitle("Mass [GeV/c^{2}]");
	92	fHistMassDimuonsOS -> Sumw2();
	93
	94	fHistMassDimuonsJpsi = new TH1D("fHistMassDimuonsJpsi", "Mass of J/#psi dimuons", 500, 0, 10);
	95	fHistMassDimuonsJpsi -> SetXTitle("Mass [GeV/c^{2}]");
	96	fHistMassDimuonsJpsi -> Sumw2();
	97
	98	fHistPtDimuonsOS = new TH1D("fHistPtDimuonsOS", "p_{T} of OS dimuons", 100, 0, 10);
	99	fHistPtDimuonsOS -> SetXTitle("p_{T} [GeV/c]");
	100	fHistPtDimuonsOS -> Sumw2();
	101
	102	fHistPtDimuonsJpsi = new TH1D("fHistPtDimuonsJpsi", "p_{T} of J/#psi dimuons", 100, 0, 10);
	103	fHistPtDimuonsJpsi -> SetXTitle("p_{T} [GeV/c]");
	104	fHistPtDimuonsJpsi -> Sumw2();
	105
	106	fHistResidualXVtxJpsi = new TH1D("fHistResidualXVtxJpsi", "J/#psi vertex residuals along x", 100, -100, 100);
	107	fHistResidualXVtxJpsi -> SetXTitle("Residual [#mum]");
	108	fHistResidualXVtxJpsi -> Sumw2();
	109
	110	fHistResidualYVtxJpsi = new TH1D("fHistResidualYVtxJpsi", "J/#psi vertex residuals along y", 100, -100, 100);
	111	fHistResidualYVtxJpsi -> SetXTitle("Residual [#mum]");
	112	fHistResidualYVtxJpsi -> Sumw2();
	113
	114	fHistResidualZVtxJpsi = new TH1D("fHistResidualZVtxJpsi", "J/#psi vertex residuals along z", 100, -1000, 1000);
	115	fHistResidualZVtxJpsi -> SetXTitle("Residual [#mum]");
	116	fHistResidualZVtxJpsi -> Sumw2();
	117
2cf4030c	118	fHistogramList->Add(fHistPtSingleMuons);
	119	fHistogramList->Add(fHistPtSingleMuonsFromJpsi);
	120	fHistogramList->Add(fHistMassDimuonsOS);
	121	fHistogramList->Add(fHistMassDimuonsJpsi);
	122	fHistogramList->Add(fHistPtDimuonsOS);
	123	fHistogramList->Add(fHistPtDimuonsJpsi);
	124	fHistogramList->Add(fHistResidualXVtxJpsi);
	125	fHistogramList->Add(fHistResidualYVtxJpsi);
	126	fHistogramList->Add(fHistResidualZVtxJpsi);
	127
	128	PostData(1, fHistogramList);
	129
d30a0b7b	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()) {
ed68e08f	212	AliAODDimuon *dimuon = new AliAODDimuon;
ed68e08f	213	dimuon->SetMuons(recMuon1,recMuon2);
bffc7f8c	214	Double_t pca[3]={0};
bffc7f8c	215	Double_t pcaQuality=0;
d30a0b7b	216	TLorentzVector kinem(0,0,0,0);
ed68e08f	217	if (!AliMFTAnalysisTools::CalculatePCA(dimuon, pca, pcaQuality, kinem)) continue;
d30a0b7b	218	fHistPtDimuonsJpsi -> Fill(kinem.Pt());
d30a0b7b	219	fHistMassDimuonsJpsi -> Fill(kinem.M());
bffc7f8c	220	fHistResidualXVtxJpsi -> Fill(1.e4*(pca[0] - fPrimaryVertex[0]));
	221	fHistResidualYVtxJpsi -> Fill(1.e4*(pca[1] - fPrimaryVertex[1]));
	222	fHistResidualZVtxJpsi -> Fill(1.e4*(pca[2] - fPrimaryVertex[2]));
ed68e08f	223	delete dimuon;
d30a0b7b	224	}
	225	}
	226	}
	227
	228	} // end of loop on 2nd muon
	229
	230	} // end of loop on 1st muon
	231
	232	//--------------------------------------------------------------------------------------------------
	233
2cf4030c	234	PostData(1, fHistogramList);
d30a0b7b	235
	236	}
	237
	238	//====================================================================================================================================================
	239
	240	void AliAnalysisTaskMFTExample::Terminate(Option_t *) {
	241
	242	// Draw result to the screen
	243	// Called once at the end of the query
	244
	245	}
	246
	247	//====================================================================================================================================================