1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
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
25 #if !defined(__CINT__) || defined(__MAKECINT__)
29 #include "AliHeader.h"
32 #include "AliFastMuonTriggerEff.h"
33 #include "AliFastMuonTrackingAcc.h"
34 #include "AliFastMuonTrackingEff.h"
35 #include "AliFastMuonTrackingRes.h"
36 #include "AliFastDetector.h"
39 #include "AliMUONTrackLight.h"
40 #include "AliMUONPairLight.h"
43 #include "TParticle.h"
46 #include "TClonesArray.h"
51 Float_t fThetaMin = 171., fThetaMax = 178.;
52 //====================================================
53 void fastMUONSim(Float_t bkg=0, char* outFileName = "fastSim_pp.root"){
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);
68 trigeff->SetCut(0); // 0... trigger pt cut at 1 GeV
69 // 1... trigger pt cut at 2 GeV
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");
77 treeOut->Branch("muons",&muonArray);
78 treeOut->Branch("dimuons",&dimuonArray);
80 // Stack of particle for each event
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");
90 RunLoader->LoadKinematics("READ");
91 Int_t nevents = RunLoader->GetNumberOfEvents();
92 printf("nevents %d\n", nevents);
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++){
99 if(iev%1000==0) printf("Event %d\n", iev);
100 RunLoader->GetEvent(iev);
101 stack = RunLoader->Stack();
102 //stack->DumpPStack();
106 muonArray->Clear(); // clean muon and dimuon arrays
107 dimuonArray->Clear();
109 Int_t nMuonsStored = 0;
110 Int_t nPart = stack->GetNtrack();
111 for(int iPart = 0; iPart < nPart; iPart++){
113 part = stack->Particle(iPart);
114 if(fabs(part->GetPdgCode()) == 13){ //muon found
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);
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();
136 //continue only if we have a muon within the MUON acceptance
137 if(theta < fThetaMin || theta > fThetaMax) continue;
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
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;
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
164 muon->FillMuonHistory(stack, part);
165 // set vertex to mother's vx,vy,vz to have the primary vertex
169 while (idMother < 0) idMother = muon->GetParentPythiaLine(id2++);
170 partMother = stack->Particle(idMother);
171 Double_t xyzMother[3] = {partMother->Vx(),
174 muon->SetVertex(xyzMother);
175 if (gRandom->Rndm() > trgeff) muon->SetTriggered(kFALSE);
176 else muon->SetTriggered(kTRUE);
177 muon->SetPxPyPz(precx,precy,precz);
179 //store the referenced track in the muonArray:
180 TClonesArray &muons = *muonArray;
181 new (muons[nMuonsStored++]) AliMUONTrackLight(*muon);
187 Int_t nmuons = muonArray->GetEntriesFast();
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);
203 RunLoader->UnloadKinematics();