]>
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 | } |