[u/mrichter/AliRoot.git] / MUON / fastMUONSim.C

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

/// \ingroup macros
/// \file fastMUONSim.C
/// \brief An example how to do the fast simulation and storing
/// the muons that survive the reconstruction
///
/// \author A. De Falco, H. Woehri, INFN Cagliari, April 2007
///
/// An example how to do the fast simulation and storing
/// the muons that survive the reconstruction in 
/// AliMUONTrackLight and AliMUONPairLight objects for
/// further analysis. Should be used together with the
/// macro fastMUONGen.C

#if !defined(__CINT__) || defined(__MAKECINT__)
//STEER includes
#include "AliStack.h"
#include "AliRun.h"
#include "AliHeader.h"

//FASTSIM includes
#include "AliFastMuonTriggerEff.h"
#include "AliFastMuonTrackingAcc.h"
#include "AliFastMuonTrackingEff.h"
#include "AliFastMuonTrackingRes.h"
#include "AliFastDetector.h"

//MUON includes
#include "AliMUONTrackLight.h"
#include "AliMUONPairLight.h"

//ROOT includes
#include "TParticle.h"
#include "TRandom.h"
#include "TFile.h"
#include "TClonesArray.h"
#include "TTree.h"

#endif

Float_t fThetaMin = 171., fThetaMax = 178.;
//====================================================
void fastMUONSim(Float_t bkg=0, char* outFileName = "fastSim_pp.root"){

  //setting up the fast simulator
  AliFastMuonTriggerEff *trigeff = new AliFastMuonTriggerEff();
  AliFastMuonTrackingAcc *acc = new AliFastMuonTrackingAcc();
  AliFastMuonTrackingEff *eff = new AliFastMuonTrackingEff();
  AliFastMuonTrackingRes *res = new AliFastMuonTrackingRes();
  acc->SetBackground(bkg);
  eff->SetBackground(bkg);
  res->SetBackground(bkg);  
  acc->Init(); 
  eff->Init(); 
  res->Init(); 
  trigeff->Init();

  trigeff->SetCut(0); // 0... trigger pt cut at 1 GeV
                      // 1... trigger pt cut at 2 GeV

  //prepare the arrays to store the generated/simulated muons
  TClonesArray *muonArray   = new TClonesArray("AliMUONTrackLight",100); 
  TClonesArray *dimuonArray = new TClonesArray("AliMUONPairLight",100); 
  TTree *treeOut = new TTree("tree","tree"); 
  TFile *fout = new TFile(outFileName,"recreate"); 
  fout->cd(); 
  treeOut->Branch("muons",&muonArray); 
  treeOut->Branch("dimuons",&dimuonArray); 
  //
  // Stack of particle for each event
  AliStack* stack;
  // Creating Run Loader and opening file containing Hits
  // AliRunLoader * RunLoader = AliRunLoader::Open("galice.root","MUONFolder","READ");
  AliRunLoader *RunLoader = AliRunLoader::Open("galice.root");
  if (RunLoader == 0x0) {
    printf(">>> Error : Error Opening %s file \n", "galice.root");
    return;
  }

  RunLoader->LoadKinematics("READ");
  Int_t nevents = RunLoader->GetNumberOfEvents();
  printf("nevents %d\n", nevents);

  TParticle *part = 0, *partMother=0;
  AliMUONTrackLight *muon = 0;
  Double_t radeg =  180./TMath::Pi(); 
  for(int iev = 0; iev < nevents; iev++){

    if(iev%1000==0) printf("Event %d\n", iev);
    RunLoader->GetEvent(iev); 
    stack = RunLoader->Stack();
    //stack->DumpPStack();

    //reset muon info
    Int_t nMuons = 0;
    muonArray->Clear();     // clean muon and dimuon arrays 
    dimuonArray->Clear(); 
    //
    Int_t nMuonsStored = 0; 
    Int_t nPart = stack->GetNtrack();
    for(int iPart = 0; iPart < nPart; iPart++){

      part = stack->Particle(iPart);
      if(fabs(part->GetPdgCode()) == 13){ //muon found

	nMuons++;
	Double_t xyz[3] = {part->Vx(),part->Vy(),part->Vz()};
	// do not take muons that were decayed after the absorber
	if(xyz[2] < -(90.+ 40.)) continue; //decay muons after 90cm + 1 int. length
	Int_t charge = (part->GetPdgCode() > 0) ? -1 : +1;
	Double_t px = part->Px();
	Double_t py = part->Py();
	Double_t pz = part->Pz();
	Double_t E = part->Energy();
	TLorentzVector pGen(px,py,pz,E);
	//
	//fast simulation code
	res->SetCharge(charge);
	eff->SetCharge(charge);
	acc->SetCharge(charge);
	Float_t p = (Float_t) pGen.P();
	Float_t pt = (Float_t) pGen.Pt();
	Float_t theta = (Float_t) radeg*pGen.Theta();
	Float_t phi = (Float_t) radeg*pGen.Phi(); 
	
	//continue only if we have a muon within the MUON acceptance
	if(theta < fThetaMin || theta > fThetaMax) continue;

	theta = 180. - theta; // correct by hand the 'new' coordinates frame
	phi   = 180. - phi; // correct by hand the 'new' coordinates frame
	Float_t prec = 0, thetarec = 0, phirec = 0; 
	res->Evaluate(p, theta ,phi, prec, thetarec, phirec);
	Float_t precx = prec * TMath::Sin(thetarec/radeg) * TMath::Cos(phirec/radeg); 
	Float_t precy = prec * TMath::Sin(thetarec/radeg) * TMath::Sin(phirec/radeg); 
	Float_t precz = prec * TMath::Cos(thetarec/radeg);
	precz = -precz; // correct by hand the 'new' coordinates frame
	precx   = -precx; // correct by hand the 'new' coordinates frame

  	Float_t trkeff = eff->Evaluate(charge, pt, theta, phi);
 	Float_t accep  = acc->Evaluate(charge, pt, theta, phi);
 	Float_t trgeff = trigeff->Evaluate(charge,pt,theta,phi);
	if(trkeff > 1.) printf("tracking efficiency > 1\n");
	if(accep > 1.) printf("acceptance efficiency > 1\n");
	if(trgeff > 1.) printf("trigger efficiency > 1\n");
	if (gRandom->Rndm() > trkeff || gRandom->Rndm() > accep) continue; 

	//only if we have a muon in the acceptance, store it as an
	//AliMUONTrackLight object and use it to form dimuons...
	muon = new AliMUONTrackLight();
	muon->SetCharge(charge);
	muon->SetTrackPDGCode(part->GetPdgCode()); //must be set, otherwise -999
	muon->SetTrackPythiaLine(iPart);//must be set, otherwise -999
	muon->SetPGen(pGen);
  	muon->FillMuonHistory(stack, part);
	// set vertex to mother's vx,vy,vz to have the primary vertex

	Int_t idMother = -1; 
	Int_t id2 = 0; 
	while (idMother < 0) idMother = muon->GetParentPythiaLine(id2++); 
	partMother = stack->Particle(idMother);
	Double_t xyzMother[3] = {partMother->Vx(),
				 partMother->Vy(),
				 partMother->Vz()};
	muon->SetVertex(xyzMother);
	if (gRandom->Rndm() > trgeff) muon->SetTriggered(kFALSE); 
	else muon->SetTriggered(kTRUE);
	muon->SetPxPyPz(precx,precy,precz);

	//store the referenced track in the muonArray:
	TClonesArray &muons = *muonArray;
	new (muons[nMuonsStored++]) AliMUONTrackLight(*muon);
	delete muon;
      }

    }//part

    Int_t nmuons = muonArray->GetEntriesFast(); 
    Int_t ndimuons = 0; 
    for(Int_t nMu1 = 0; nMu1 < nmuons-1; nMu1++) { 
      AliMUONTrackLight* mu1 = (AliMUONTrackLight*) muonArray->At(nMu1); 
      for(Int_t nMu2 = nMu1+1; nMu2 < nmuons; nMu2++){
	AliMUONTrackLight* mu2 = (AliMUONTrackLight*) muonArray->At(nMu2); 
	AliMUONPairLight *dimuLight = new AliMUONPairLight(); 
	dimuLight->SetMuons(*mu1, *mu2);
	TClonesArray &dimuons = *dimuonArray;
	new (dimuons[ndimuons++]) AliMUONPairLight(*dimuLight);
      }
    }// end dimuons
    treeOut->Fill(); 
    stack->Reset();
  }//end of events

  RunLoader->UnloadKinematics();
  fout->cd(); 
  treeOut->Write();
}
Commit	Line	Data
9c090011	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
e54bf126	16	/* $Id$ */
9c090011	17
e54bf126	18	/// \ingroup macros
	19	/// \file fastMUONSim.C
	20	/// \brief An example how to do the fast simulation and storing
	21	/// the muons that survive the reconstruction
	22	///
	23	/// \author A. De Falco, H. Woehri, INFN Cagliari, April 2007
	24	///
	25	/// An example how to do the fast simulation and storing
	26	/// the muons that survive the reconstruction in
	27	/// AliMUONTrackLight and AliMUONPairLight objects for
	28	/// further analysis. Should be used together with the
	29	/// macro fastMUONGen.C
9c090011	30
	31	#if !defined(__CINT__) \|\| defined(__MAKECINT__)
	32	//STEER includes
	33	#include "AliStack.h"
	34	#include "AliRun.h"
	35	#include "AliHeader.h"
	36
	37	//FASTSIM includes
	38	#include "AliFastMuonTriggerEff.h"
	39	#include "AliFastMuonTrackingAcc.h"
	40	#include "AliFastMuonTrackingEff.h"
	41	#include "AliFastMuonTrackingRes.h"
	42	#include "AliFastDetector.h"
	43
	44	//MUON includes
	45	#include "AliMUONTrackLight.h"
	46	#include "AliMUONPairLight.h"
	47
	48	//ROOT includes
	49	#include "TParticle.h"
	50	#include "TRandom.h"
	51	#include "TFile.h"
	52	#include "TClonesArray.h"
	53	#include "TTree.h"
	54
	55	#endif
	56
	57	Float_t fThetaMin = 171., fThetaMax = 178.;
	58	//====================================================
	59	void fastMUONSim(Float_t bkg=0, char* outFileName = "fastSim_pp.root"){
	60
	61	//setting up the fast simulator
	62	AliFastMuonTriggerEff *trigeff = new AliFastMuonTriggerEff();
	63	AliFastMuonTrackingAcc *acc = new AliFastMuonTrackingAcc();
	64	AliFastMuonTrackingEff *eff = new AliFastMuonTrackingEff();
	65	AliFastMuonTrackingRes *res = new AliFastMuonTrackingRes();
	66	acc->SetBackground(bkg);
	67	eff->SetBackground(bkg);
	68	res->SetBackground(bkg);
	69	acc->Init();
	70	eff->Init();
	71	res->Init();
	72	trigeff->Init();
	73
	74	trigeff->SetCut(0); // 0... trigger pt cut at 1 GeV
	75	// 1... trigger pt cut at 2 GeV
	76
	77	//prepare the arrays to store the generated/simulated muons
	78	TClonesArray *muonArray = new TClonesArray("AliMUONTrackLight",100);
	79	TClonesArray *dimuonArray = new TClonesArray("AliMUONPairLight",100);
	80	TTree *treeOut = new TTree("tree","tree");
	81	TFile *fout = new TFile(outFileName,"recreate");
	82	fout->cd();
	83	treeOut->Branch("muons",&muonArray);
	84	treeOut->Branch("dimuons",&dimuonArray);
	85	//
	86	// Stack of particle for each event
	87	AliStack* stack;
	88	// Creating Run Loader and opening file containing Hits
	89	// AliRunLoader * RunLoader = AliRunLoader::Open("galice.root","MUONFolder","READ");
	90	AliRunLoader *RunLoader = AliRunLoader::Open("galice.root");
	91	if (RunLoader == 0x0) {
	92	printf(">>> Error : Error Opening %s file \n", "galice.root");
	93	return;
94	}
95
96	RunLoader->LoadKinematics("READ");
97	Int_t nevents = RunLoader->GetNumberOfEvents();
98	printf("nevents %d\n", nevents);
99
100	TParticle part = 0, partMother=0;
101	AliMUONTrackLight *muon = 0;
102	Double_t radeg = 180./TMath::Pi();
103	for(int iev = 0; iev < nevents; iev++){
104
105	if(iev%1000==0) printf("Event %d\n", iev);
106	RunLoader->GetEvent(iev);
107	stack = RunLoader->Stack();
108	//stack->DumpPStack();
109
110	//reset muon info
111	Int_t nMuons = 0;
112	muonArray->Clear(); // clean muon and dimuon arrays
113	dimuonArray->Clear();
114	//
115	Int_t nMuonsStored = 0;
116	Int_t nPart = stack->GetNtrack();
117	for(int iPart = 0; iPart < nPart; iPart++){
118
119	part = stack->Particle(iPart);
120	if(fabs(part->GetPdgCode()) == 13){ //muon found
121
122	nMuons++;
123	Double_t xyz[3] = {part->Vx(),part->Vy(),part->Vz()};
124	// do not take muons that were decayed after the absorber
125	if(xyz[2] < -(90.+ 40.)) continue; //decay muons after 90cm + 1 int. length
126	Int_t charge = (part->GetPdgCode() > 0) ? -1 : +1;
127	Double_t px = part->Px();
128	Double_t py = part->Py();
129	Double_t pz = part->Pz();
130	Double_t E = part->Energy();
131	TLorentzVector pGen(px,py,pz,E);
132	//
133	//fast simulation code
134	res->SetCharge(charge);
135	eff->SetCharge(charge);
136	acc->SetCharge(charge);
137	Float_t p = (Float_t) pGen.P();
138	Float_t pt = (Float_t) pGen.Pt();
139	Float_t theta = (Float_t) radeg*pGen.Theta();
140	Float_t phi = (Float_t) radeg*pGen.Phi();
141
142	//continue only if we have a muon within the MUON acceptance
143	if(theta < fThetaMin \|\| theta > fThetaMax) continue;
144
145	theta = 180. - theta; // correct by hand the 'new' coordinates frame
146	phi = 180. - phi; // correct by hand the 'new' coordinates frame
147	Float_t prec = 0, thetarec = 0, phirec = 0;
148	res->Evaluate(p, theta ,phi, prec, thetarec, phirec);
149	Float_t precx = prec * TMath::Sin(thetarec/radeg) * TMath::Cos(phirec/radeg);
150	Float_t precy = prec * TMath::Sin(thetarec/radeg) * TMath::Sin(phirec/radeg);
151	Float_t precz = prec * TMath::Cos(thetarec/radeg);
152	precz = -precz; // correct by hand the 'new' coordinates frame
153	precx = -precx; // correct by hand the 'new' coordinates frame
154
155	Float_t trkeff = eff->Evaluate(charge, pt, theta, phi);
156	Float_t accep = acc->Evaluate(charge, pt, theta, phi);
157	Float_t trgeff = trigeff->Evaluate(charge,pt,theta,phi);
158	if(trkeff > 1.) printf("tracking efficiency > 1\n");
159	if(accep > 1.) printf("acceptance efficiency > 1\n");
160	if(trgeff > 1.) printf("trigger efficiency > 1\n");
161	if (gRandom->Rndm() > trkeff \|\| gRandom->Rndm() > accep) continue;
162
163	//only if we have a muon in the acceptance, store it as an
164	//AliMUONTrackLight object and use it to form dimuons...
165	muon = new AliMUONTrackLight();
166	muon->SetCharge(charge);
167	muon->SetTrackPDGCode(part->GetPdgCode()); //must be set, otherwise -999
168	muon->SetTrackPythiaLine(iPart);//must be set, otherwise -999
169	muon->SetPGen(pGen);
170	muon->FillMuonHistory(stack, part);
171	// set vertex to mother's vx,vy,vz to have the primary vertex
172
173	Int_t idMother = -1;
174	Int_t id2 = 0;
175	while (idMother < 0) idMother = muon->GetParentPythiaLine(id2++);
176	partMother = stack->Particle(idMother);
177	Double_t xyzMother[3] = {partMother->Vx(),
178	partMother->Vy(),
179	partMother->Vz()};
180	muon->SetVertex(xyzMother);
181	if (gRandom->Rndm() > trgeff) muon->SetTriggered(kFALSE);
182	else muon->SetTriggered(kTRUE);
183	muon->SetPxPyPz(precx,precy,precz);
184
185	//store the referenced track in the muonArray:
186	TClonesArray &muons = *muonArray;
187	new (muons[nMuonsStored++]) AliMUONTrackLight(*muon);
188	delete muon;
189	}
190
191	}//part
192
193	Int_t nmuons = muonArray->GetEntriesFast();
194	Int_t ndimuons = 0;
195	for(Int_t nMu1 = 0; nMu1 < nmuons-1; nMu1++) {
196	AliMUONTrackLight* mu1 = (AliMUONTrackLight*) muonArray->At(nMu1);
197	for(Int_t nMu2 = nMu1+1; nMu2 < nmuons; nMu2++){
198	AliMUONTrackLight* mu2 = (AliMUONTrackLight*) muonArray->At(nMu2);
199	AliMUONPairLight *dimuLight = new AliMUONPairLight();
200	dimuLight->SetMuons(mu1, mu2);
201	TClonesArray &dimuons = *dimuonArray;
202	new (dimuons[ndimuons++]) AliMUONPairLight(*dimuLight);
203	}
204	}// end dimuons
205	treeOut->Fill();
206	stack->Reset();
207	}//end of events
208
209	RunLoader->UnloadKinematics();
210	fout->cd();
211	treeOut->Write();
212	}