Separating run-dependent mapping data from data, which are not
[u/mrichter/AliRoot.git] / MUON / fastMUONSim.C
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 /// \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
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 }