[u/mrichter/AliRoot.git] / PWG3 / muon / AliAnalysisTaskMuonAODfromGeneral.cxx

#define AliAnalysisTaskMuonAODfromGeneral_cxx

// 19 Nov 2007
// Class implementation for the specific muon AOD generation
// Extracts only muon tracks from a general AOD and builds dimuons
// Livio Bianchi, Universita' di Torino


#include "TTree.h"
#include "TROOT.h"
#include "TH1.h"
#include "TCanvas.h"
#include "TSystem.h"
#include "TRandom.h"

#include "AliAODEvent.h"
#include "AliAnalysisTask.h"
#include "AliAnalysisTaskMuonAODfromGeneral.h"
#include "AliAODHandler.h"
#include "AliAnalysisManager.h"
#include "AliAODEventInfo.h"
#include "AliAODDimuon.h"

ClassImp(AliAnalysisTaskMuonAODfromGeneral)

//________________________________________________________________________
AliAnalysisTaskMuonAODfromGeneral::AliAnalysisTaskMuonAODfromGeneral(const char *name, Double_t BeamEnergy):
	AliAnalysisTask(name, "AnalysisTaskMuonAODfromGeneral"),
	fInfos(0),
	fDimuons(0),
        fChain(0),
        fOrgAOD(0), 
	fNewAOD(0), 
	ft(0),
	fBeamEnergy(0)
{  // Constructor.
  // Input slot #0 works with an Ntuple
  SetBeamEnergy(BeamEnergy);
  DefineInput(0, TChain::Class());
  // Output slot #0 writes into a TTree container
  DefineOutput(0, TTree::Class());
}

//________________________________________________________________________
void AliAnalysisTaskMuonAODfromGeneral::ConnectInputData(Option_t *) {
  printf("   ConnectInputData %s\n", GetName());
  // New connection to the AOD
  fChain = (TChain*)GetInputData(0);
  fOrgAOD = new AliAODEvent();
  fOrgAOD->ReadFromTree(fChain); // This checks also if the branch address is already set
}

//________________________________________________________________________
void AliAnalysisTaskMuonAODfromGeneral::CreateOutputObjects() {
  printf("Creating output objects\n");
  if(!fNewAOD)fNewAOD=new AliAODEvent();
  fNewAOD->CreateStdContent();
 
  if (!ft) {
    ft = new TTree("AOD","Muon AOD tree",0);
    ft->BranchRef();
    ft->Branch(fNewAOD->GetList());
    //ft->Print();
 }
 
  //Add new Objects: MuonInfos & Dimuons
   fInfos = new AliAODEventInfo();
   fNewAOD->AddObject(fInfos);
   const char *nameInfos = "MuonInfos";
   fInfos->SetName(nameInfos);
   ft->Branch(nameInfos, &fInfos);
   fDimuons = new TClonesArray("AliAODDimuon",0);
   fNewAOD->AddObject(fDimuons);
   const char *nameDimuons = "Dimuons";
   fDimuons->SetName(nameDimuons);
   ft->Branch(nameDimuons, &fDimuons);
      
}

//________________________________________________________________________
void AliAnalysisTaskMuonAODfromGeneral::Exec(Option_t *) {
  static int ncall=0;
  // Task making a Muon AOD
  // Get input data
  TChain *chain = (TChain*)GetInputData(0);
  Long64_t ientry = chain->GetReadEntry();
  ientry=ientry;

  if (!fOrgAOD) return;
  
  Int_t nTracks=fOrgAOD->GetNumberOfTracks();
  Int_t nMuTracks=0;
  Int_t nPosTracks = 0;
  Int_t mutrNumb[10]; for (Int_t i=0; i<10; mutrNumb[i++]=0);
  for (Int_t iTrack=0; iTrack<nTracks; iTrack++){
       const Double_t *trackpid=fOrgAOD->GetTrack(iTrack)->PID();
       if (trackpid[AliAODTrack::kMuon]==1.) {
       		mutrNumb[nMuTracks]=iTrack;
		nMuTracks++;
		}
       if (fOrgAOD->GetTrack(iTrack)->Charge()> 0) nPosTracks++;
       }
  
  Bool_t ExistMuon=kFALSE;
  if (nMuTracks>0) ExistMuon=kTRUE;
  
  //--------------------------------------------------------------------
  // set arrays and pointers
  Double_t pos[3];
  Double_t covVtx[6];
  
    // Access to the header
    AliAODHeader *header = fNewAOD->GetHeader();

  // fill the header
  header->SetRunNumber       (fOrgAOD->GetRunNumber()       );
  header->SetBunchCrossNumber(fOrgAOD->GetBunchCrossNumber());
  header->SetOrbitNumber     (fOrgAOD->GetOrbitNumber()     );
  header->SetPeriodNumber    (fOrgAOD->GetPeriodNumber()    );
  header->SetTriggerMask     (fOrgAOD->GetTriggerMask()     ); 
  header->SetTriggerCluster  (fOrgAOD->GetTriggerCluster()  );
  header->SetEventType       (fOrgAOD->GetEventType()       );
  header->SetMagneticField   (fOrgAOD->GetMagneticField()   );
  header->SetZDCN1Energy     (fOrgAOD->GetZDCN1Energy()     );
  header->SetZDCP1Energy     (fOrgAOD->GetZDCP1Energy()     );
  header->SetZDCN2Energy     (fOrgAOD->GetZDCN2Energy()     );
  header->SetZDCP2Energy     (fOrgAOD->GetZDCP2Energy()     );
  header->SetZDCEMEnergy     (0,fOrgAOD->GetZDCEMEnergy(0)  );
  header->SetZDCEMEnergy     (1,fOrgAOD->GetZDCEMEnergy(1)  );
  header->SetRefMultiplicity   (nTracks);
  header->SetRefMultiplicityPos(nPosTracks);
  header->SetRefMultiplicityNeg(nTracks - nPosTracks);
  header->SetMuonMagFieldScale(-999.); // FIXME
  header->SetCentrality(-999.);        // FIXME

  Int_t nVertices = 0;
  const AliAODVertex *vtx = fOrgAOD->GetPrimaryVertex();
  if(vtx)nVertices = 1;

  Int_t nDimuons=nMuTracks*(nMuTracks-1)/2;
  
  // Access to the AOD container of vertices
  TClonesArray &vertices = *(fNewAOD->GetVertices());
  Int_t jVertices=0;

  // Access to the AOD container of tracks
   TClonesArray &tracks = *(fNewAOD->GetTracks());
   tracks.Delete();
   Int_t jTracks=0; 


  vtx->GetCovMatrix(covVtx); //covariance matrix

  AliAODVertex * primary = new(vertices[jVertices++])
      AliAODVertex(pos, covVtx, vtx->GetChi2perNDF(), NULL, AliAODVertex::kPrimary);

  static int ncal=0;
  static int numtracks=0;
  for (Int_t iTrack=0; iTrack<nMuTracks; iTrack++) {
	AliAODTrack *aodTrack;
	primary->AddDaughter( aodTrack=new(tracks[jTracks++])
	AliAODTrack((*(fOrgAOD->GetTrack(mutrNumb[iTrack])))));
        ncal++;
  }
  numtracks+=jTracks;

  fDimuons->Delete();
  fDimuons->Expand(nDimuons);
  TClonesArray &rDimuons = *fDimuons;
    
  Int_t jDimuons=0;
  for(Int_t i=0; i<nMuTracks; i++){
     for(Int_t j=i+1; j<nMuTracks; j++){
       new(rDimuons[jDimuons++]) AliAODDimuon(tracks[i],tracks[j],fInfos);
     } 
  }
  
  fInfos->SetBeamEnergy(fBeamEnergy);
  fInfos->SetEv(fNewAOD);
  fInfos->SetEi(fInfos);
  fInfos->SetHe(header);
  fInfos->SetTr(fNewAOD->GetTracks());
  fInfos->SetDi(fDimuons);
  fInfos->SelectTriggerBits(0,1,2,3,4,5);
  if(ExistMuon) ft->Fill();
  ncall++;

  // Post final data. It will be written to a file with option "RECREATE"
  PostData(0, ft);
}      

//________________________________________________________________________
void AliAnalysisTaskMuonAODfromGeneral::Terminate(Option_t *) {
  // Draw some histogram at the end.
  ft->GetUserInfo()->Add(fNewAOD);
  ft->Write("",TObject::kOverwrite);
  ft->ResetBranchAddresses();
  if (!gROOT->IsBatch()) {
    TCanvas *c1 = new TCanvas("c1","Pt",10,10,310,310);
    c1->SetFillColor(10);
    c1->SetHighLightColor(10);
    
    c1->cd(1)->SetLeftMargin(0.15);
    c1->cd(1)->SetBottomMargin(0.15);  
    c1->cd(1)->SetLogy();
    ft->Draw("tracks.Pt()");
    //ft->Draw("dimuons.M()");
  }
}
Commit	Line	Data
71b7d225	1	#define AliAnalysisTaskMuonAODfromGeneral_cxx
	2
	3	// 19 Nov 2007
aba89748	4	// Class implementation for the specific muon AOD generation
	5	// Extracts only muon tracks from a general AOD and builds dimuons
	6	// Livio Bianchi, Universita' di Torino
	7
71b7d225	8
	9	#include "TTree.h"
	10	#include "TROOT.h"
	11	#include "TH1.h"
	12	#include "TCanvas.h"
	13	#include "TSystem.h"
	14	#include "TRandom.h"
	15
	16	#include "AliAODEvent.h"
aba89748	17	#include "AliAnalysisTask.h"
71b7d225	18	#include "AliAnalysisTaskMuonAODfromGeneral.h"
	19	#include "AliAODHandler.h"
	20	#include "AliAnalysisManager.h"
	21	#include "AliAODEventInfo.h"
	22	#include "AliAODDimuon.h"
	23
	24	ClassImp(AliAnalysisTaskMuonAODfromGeneral)
	25
	26	//________________________________________________________________________
	27	AliAnalysisTaskMuonAODfromGeneral::AliAnalysisTaskMuonAODfromGeneral(const char *name, Double_t BeamEnergy):
	28	AliAnalysisTask(name, "AnalysisTaskMuonAODfromGeneral"),
	29	fInfos(0),
	30	fDimuons(0),
8b8b0b7a	31	fChain(0),
8b8b0b7a	32	fOrgAOD(0),
71b7d225	33	fNewAOD(0),
	34	ft(0),
	35	fBeamEnergy(0)
	36	{ // Constructor.
	37	// Input slot #0 works with an Ntuple
	38	SetBeamEnergy(BeamEnergy);
	39	DefineInput(0, TChain::Class());
	40	// Output slot #0 writes into a TTree container
	41	DefineOutput(0, TTree::Class());
	42	}
	43
	44	//________________________________________________________________________
	45	void AliAnalysisTaskMuonAODfromGeneral::ConnectInputData(Option_t *) {
	46	printf(" ConnectInputData %s\n", GetName());
	47	// New connection to the AOD
	48	fChain = (TChain*)GetInputData(0);
	49	fOrgAOD = new AliAODEvent();
	50	fOrgAOD->ReadFromTree(fChain); // This checks also if the branch address is already set
	51	}
	52
	53	//________________________________________________________________________
	54	void AliAnalysisTaskMuonAODfromGeneral::CreateOutputObjects() {
	55	printf("Creating output objects\n");
	56	if(!fNewAOD)fNewAOD=new AliAODEvent();
	57	fNewAOD->CreateStdContent();
	58
	59	if (!ft) {
	60	ft = new TTree("AOD","Muon AOD tree",0);
	61	ft->BranchRef();
	62	ft->Branch(fNewAOD->GetList());
	63	//ft->Print();
	64	}
	65
	66	//Add new Objects: MuonInfos & Dimuons
	67	fInfos = new AliAODEventInfo();
	68	fNewAOD->AddObject(fInfos);
	69	const char *nameInfos = "MuonInfos";
	70	fInfos->SetName(nameInfos);
	71	ft->Branch(nameInfos, &fInfos);
	72	fDimuons = new TClonesArray("AliAODDimuon",0);
	73	fNewAOD->AddObject(fDimuons);
	74	const char *nameDimuons = "Dimuons";
	75	fDimuons->SetName(nameDimuons);
	76	ft->Branch(nameDimuons, &fDimuons);
	77
	78	}
	79
	80	//________________________________________________________________________
	81	void AliAnalysisTaskMuonAODfromGeneral::Exec(Option_t *) {
	82	static int ncall=0;
	83	// Task making a Muon AOD
	84	// Get input data
	85	TChain chain = (TChain)GetInputData(0);
	86	Long64_t ientry = chain->GetReadEntry();
	87	ientry=ientry;
	88
	89	if (!fOrgAOD) return;
	90
	91	Int_t nTracks=fOrgAOD->GetNumberOfTracks();
	92	Int_t nMuTracks=0;
	93	Int_t nPosTracks = 0;
	94	Int_t mutrNumb[10]; for (Int_t i=0; i<10; mutrNumb[i++]=0);
	95	for (Int_t iTrack=0; iTrack<nTracks; iTrack++){
	96	const Double_t *trackpid=fOrgAOD->GetTrack(iTrack)->PID();
97	if (trackpid[AliAODTrack::kMuon]==1.) {
98	mutrNumb[nMuTracks]=iTrack;
99	nMuTracks++;
100	}
101	if (fOrgAOD->GetTrack(iTrack)->Charge()> 0) nPosTracks++;
102	}
103
104	Bool_t ExistMuon=kFALSE;
105	if (nMuTracks>0) ExistMuon=kTRUE;
106
107	//--------------------------------------------------------------------
108	// set arrays and pointers
109	Double_t pos[3];
110	Double_t covVtx[6];
111
112	// Access to the header
113	AliAODHeader *header = fNewAOD->GetHeader();
114
115	// fill the header
116	header->SetRunNumber (fOrgAOD->GetRunNumber() );
117	header->SetBunchCrossNumber(fOrgAOD->GetBunchCrossNumber());
118	header->SetOrbitNumber (fOrgAOD->GetOrbitNumber() );
119	header->SetPeriodNumber (fOrgAOD->GetPeriodNumber() );
120	header->SetTriggerMask (fOrgAOD->GetTriggerMask() );
121	header->SetTriggerCluster (fOrgAOD->GetTriggerCluster() );
122	header->SetEventType (fOrgAOD->GetEventType() );
123	header->SetMagneticField (fOrgAOD->GetMagneticField() );
124	header->SetZDCN1Energy (fOrgAOD->GetZDCN1Energy() );
125	header->SetZDCP1Energy (fOrgAOD->GetZDCP1Energy() );
126	header->SetZDCN2Energy (fOrgAOD->GetZDCN2Energy() );
127	header->SetZDCP2Energy (fOrgAOD->GetZDCP2Energy() );
128	header->SetZDCEMEnergy (0,fOrgAOD->GetZDCEMEnergy(0) );
129	header->SetZDCEMEnergy (1,fOrgAOD->GetZDCEMEnergy(1) );
130	header->SetRefMultiplicity (nTracks);
131	header->SetRefMultiplicityPos(nPosTracks);
132	header->SetRefMultiplicityNeg(nTracks - nPosTracks);
133	header->SetMuonMagFieldScale(-999.); // FIXME
134	header->SetCentrality(-999.); // FIXME
135
136	Int_t nVertices = 0;
137	const AliAODVertex *vtx = fOrgAOD->GetPrimaryVertex();
138	if(vtx)nVertices = 1;
139
140	Int_t nDimuons=nMuTracks*(nMuTracks-1)/2;
141
142	// Access to the AOD container of vertices
143	TClonesArray &vertices = *(fNewAOD->GetVertices());
144	Int_t jVertices=0;
145
146	// Access to the AOD container of tracks
147	TClonesArray &tracks = *(fNewAOD->GetTracks());
148	tracks.Delete();
149	Int_t jTracks=0;
150
151
152	vtx->GetCovMatrix(covVtx); //covariance matrix
153
154	AliAODVertex * primary = new(vertices[jVertices++])
155	AliAODVertex(pos, covVtx, vtx->GetChi2perNDF(), NULL, AliAODVertex::kPrimary);
156
157	static int ncal=0;
158	static int numtracks=0;
159	for (Int_t iTrack=0; iTrack<nMuTracks; iTrack++) {
160	AliAODTrack *aodTrack;
161	primary->AddDaughter( aodTrack=new(tracks[jTracks++])
162	AliAODTrack((*(fOrgAOD->GetTrack(mutrNumb[iTrack])))));
163	ncal++;
164	}
165	numtracks+=jTracks;
166
167	fDimuons->Delete();
168	fDimuons->Expand(nDimuons);
169	TClonesArray &rDimuons = *fDimuons;
170
171	Int_t jDimuons=0;
172	for(Int_t i=0; i<nMuTracks; i++){
173	for(Int_t j=i+1; j<nMuTracks; j++){
174	new(rDimuons[jDimuons++]) AliAODDimuon(tracks[i],tracks[j],fInfos);
175	}
176	}
177
178	fInfos->SetBeamEnergy(fBeamEnergy);
aba89748	179	fInfos->SetEv(fNewAOD);
	180	fInfos->SetEi(fInfos);
	181	fInfos->SetHe(header);
	182	fInfos->SetTr(fNewAOD->GetTracks());
	183	fInfos->SetDi(fDimuons);
	184	fInfos->SelectTriggerBits(0,1,2,3,4,5);
71b7d225	185	if(ExistMuon) ft->Fill();
	186	ncall++;
	187
	188	// Post final data. It will be written to a file with option "RECREATE"
	189	PostData(0, ft);
	190	}
	191
	192	//________________________________________________________________________
	193	void AliAnalysisTaskMuonAODfromGeneral::Terminate(Option_t *) {
	194	// Draw some histogram at the end.
	195	ft->GetUserInfo()->Add(fNewAOD);
	196	ft->Write("",TObject::kOverwrite);
	197	ft->ResetBranchAddresses();
	198	if (!gROOT->IsBatch()) {
	199	TCanvas *c1 = new TCanvas("c1","Pt",10,10,310,310);
	200	c1->SetFillColor(10);
	201	c1->SetHighLightColor(10);
	202
	203	c1->cd(1)->SetLeftMargin(0.15);
	204	c1->cd(1)->SetBottomMargin(0.15);
	205	c1->cd(1)->SetLogy();
	206	ft->Draw("tracks.Pt()");
	207	//ft->Draw("dimuons.M()");
	208	}
	209	}