Fixes for #86059: Install data when ALICE_ROOT does not point to source (Christian)
[u/mrichter/AliRoot.git] / MUON / fastMUONSim.C
CommitLineData
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
57Float_t fThetaMin = 171., fThetaMax = 178.;
58//====================================================
59void 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}