[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: */

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