New macros to run fast simulation
[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
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
51Float_t fThetaMin = 171., fThetaMax = 178.;
52//====================================================
53void 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}