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