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 **************************************************************************/
19 /// \file fastMUONSim.C
20 /// \brief An example how to do the fast simulation and storing
21 /// the muons that survive the reconstruction
23 /// \author A. De Falco, H. Woehri, INFN Cagliari, April 2007
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
31 #if !defined(__CINT__) || defined(__MAKECINT__)
35 #include "AliHeader.h"
38 #include "AliFastMuonTriggerEff.h"
39 #include "AliFastMuonTrackingAcc.h"
40 #include "AliFastMuonTrackingEff.h"
41 #include "AliFastMuonTrackingRes.h"
42 #include "AliFastDetector.h"
45 #include "AliMUONTrackLight.h"
46 #include "AliMUONPairLight.h"
49 #include "TParticle.h"
52 #include "TClonesArray.h"
57 Float_t fThetaMin = 171., fThetaMax = 178.;
58 //====================================================
59 void fastMUONSim(Float_t bkg=0, char* outFileName = "fastSim_pp.root"){
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);
74 trigeff->SetCut(0); // 0... trigger pt cut at 1 GeV
75 // 1... trigger pt cut at 2 GeV
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");
83 treeOut->Branch("muons",&muonArray);
84 treeOut->Branch("dimuons",&dimuonArray);
86 // Stack of particle for each event
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");
96 RunLoader->LoadKinematics("READ");
97 Int_t nevents = RunLoader->GetNumberOfEvents();
98 printf("nevents %d\n", nevents);
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++){
105 if(iev%1000==0) printf("Event %d\n", iev);
106 RunLoader->GetEvent(iev);
107 stack = RunLoader->Stack();
108 //stack->DumpPStack();
112 muonArray->Clear(); // clean muon and dimuon arrays
113 dimuonArray->Clear();
115 Int_t nMuonsStored = 0;
116 Int_t nPart = stack->GetNtrack();
117 for(int iPart = 0; iPart < nPart; iPart++){
119 part = stack->Particle(iPart);
120 if(fabs(part->GetPdgCode()) == 13){ //muon found
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);
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();
142 //continue only if we have a muon within the MUON acceptance
143 if(theta < fThetaMin || theta > fThetaMax) continue;
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
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;
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
170 muon->FillMuonHistory(stack, part);
171 // set vertex to mother's vx,vy,vz to have the primary vertex
175 while (idMother < 0) idMother = muon->GetParentPythiaLine(id2++);
176 partMother = stack->Particle(idMother);
177 Double_t xyzMother[3] = {partMother->Vx(),
180 muon->SetVertex(xyzMother);
181 if (gRandom->Rndm() > trgeff) muon->SetTriggered(kFALSE);
182 else muon->SetTriggered(kTRUE);
183 muon->SetPxPyPz(precx,precy,precz);
185 //store the referenced track in the muonArray:
186 TClonesArray &muons = *muonArray;
187 new (muons[nMuonsStored++]) AliMUONTrackLight(*muon);
193 Int_t nmuons = muonArray->GetEntriesFast();
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);
209 RunLoader->UnloadKinematics();