Reading muon trigger scalers with the DA of the muon trigger and transfer
[u/mrichter/AliRoot.git] / FASTSIM / fastGenQuarkonia.C
1 AliGenerator*  CreateGenerator();
2
3 void fastGenQuarkonia(Int_t nev = 1, char* filename = "galice.root")
4 {
5 //  Update data base
6     AliPDG::AddParticlesToPdgDataBase();
7     
8 // Create the fast tracker
9     AliFastMuonTriggerEff *trigeff = new AliFastMuonTriggerEff();
10     AliFastMuonTrackingAcc *acc = new AliFastMuonTrackingAcc();
11     AliFastMuonTrackingEff *eff = new AliFastMuonTrackingEff();
12     AliFastMuonTrackingRes *res = new AliFastMuonTrackingRes();
13     acc->SetBackground(0);
14     eff->SetBackground(0);
15     res->SetBackground(0);  
16     acc->Init(); 
17     eff->Init(); 
18     res->Init(); 
19     AliFastDetector* tracker = new AliFastDetector("Tracker", "Muon Tracker");
20     tracker->AddResponse(trigeff);
21     tracker->AddResponse(acc);
22     tracker->AddResponse(eff);
23     tracker->AddResponse(res);
24     tracker->Init();
25     tracker->Dump();
26 //  Histos
27     TH1F* massHU = new TH1F("massHU", "Mass Spectrum:                ", 500, 5, 15.);
28     TH1F* massHS = new TH1F("massHS", "Mass Spectrum Smeared:        ", 500, 5, 15.);
29     TH2F* etaH   = new TH2F("etaH",   "eta vs etas        ", 100, 0., 5., 100, 0., 5.);
30     TH1F* etadH  = new TH1F("etaH",   "eta vs etas        ", 100, -1., 1.);
31 //
32 //                        Construction
33 //
34     AliRunLoader* rl = AliRunLoader::Open("galice.root","FASTRUN","recreate");
35     
36     rl->SetCompressionLevel(2);
37     rl->SetNumberOfEventsPerFile(nev);
38     rl->LoadKinematics("RECREATE");
39     rl->MakeTree("E");
40     gAlice->SetRunLoader(rl);
41
42 //  Create stack
43     rl->MakeStack();
44     AliStack* stack      = rl->Stack();
45  
46 //  Header
47     AliHeader* header = rl->GetHeader();
48 //
49 //  Create and Initialize Generator
50     AliGenerator *gener = CreateGenerator();
51     gener->Init();
52     gener->SetStack(stack);
53     
54 //
55 //                        Event Loop
56 //
57     Int_t iev;
58      
59     for (iev = 0; iev < nev; iev++) {
60
61         printf("\n \n Event number %d \n \n", iev);
62         
63
64 //  Initialize event
65         header->Reset(0,iev);
66         rl->SetEventNumber(iev);
67         stack->Reset();
68         rl->MakeTree("K");
69 //  Generate event
70         gener->Generate();
71                 
72 //  Analysis
73         Int_t npart = stack->GetNprimary();
74         for (Int_t part=0; part<npart; part++) {
75             TParticle *MPart = stack->Particle(part);
76             Int_t mpart  = MPart->GetPdgCode();
77             if (mpart != 553 && mpart != 100553 && mpart != 200553) continue;
78             Int_t ch1 = MPart->GetFirstDaughter();
79             Int_t ch2 = MPart->GetLastDaughter();
80             
81             TParticle *muon1 = stack->Particle(ch1);
82             TParticle *muon2 = stack->Particle(ch2);
83             
84             Float_t theta1 = muon1->Theta();
85             Float_t thetad1 = theta1 * 180./TMath::Pi();
86             Float_t eta1   = muon1->Eta();
87             Float_t pt1    = muon1->Pt();
88             Float_t phi1   = muon1->Phi();
89             Float_t phid1  = phi1 * 180./TMath::Pi() - 180.;
90             Float_t p1     = muon1->P();
91
92             if (thetad1 > 9. || thetad1 < 2.) continue;
93
94             Float_t theta2 = muon2->Theta();
95             Float_t thetad2 = theta2 * 180./TMath::Pi();
96             Float_t eta2   = muon2->Eta();
97             Float_t pt2    = muon2->Pt();
98             Float_t phi2   = muon2->Phi();
99             Float_t phid2  = phi2 * 180./TMath::Pi() - 180.;
100             Float_t p2     = muon2->P();
101
102
103             if (thetad2 > 9. || thetad2 < 2.) continue;     
104             
105             Float_t dphi   = phi1 - phi2;
106             Float_t deta   = eta1 - eta2;
107             Float_t m      =  TMath::Sqrt(2. * pt1 * pt2 * (TMath::CosH(deta) - TMath::Cos(dphi)));
108
109             massHU->Fill(m);
110 // Smear
111             // the mu+
112             Float_t thetas1, phis1, ps1, thetas2, phis2, ps2, pts1, pts2, etas1, etas2;
113             Float_t trigeffpL, trigeffpH, trigeffnL, trigeffnH, trigeffnA;
114             
115             res->SetCharge(1);
116             eff->SetCharge(1);
117             acc->SetCharge(1);
118             res->Evaluate(p1, thetad1, phid1, ps1, thetas1, phis1);
119
120             
121             Float_t effp     = eff->Evaluate(pt1, thetad1, phid1);
122             Float_t accp     = acc->Evaluate(pt1, thetad1, phid1);
123             trigeff->Evaluate(1, pt1, thetad1, phid1, trigeffpL, trigeffpH, trigeffnA);
124             thetas1 *= TMath::Pi()/180.;
125             phis1 *= TMath::Pi()/180.;
126             Float_t etas;
127             
128             if (TMath::Tan(thetas1/2.) > 0.) {
129                 etas  = -TMath::Log(TMath::Tan(thetas1/2.));
130                 etaH->Fill(etas, eta1);
131                 etadH->Fill(etas-eta1);
132             }
133             
134             // the mu- 
135             res->SetCharge(-1);
136             acc->SetCharge(-1);
137             eff->SetCharge(-1);
138             res->Evaluate(p2, thetad2, phid2, ps2, thetas2, phis2);
139             Float_t effn     = eff->Evaluate(pt2, thetad2, phid2);
140             Float_t accn     = acc->Evaluate(pt2, thetad2, phid2);
141             trigeff->Evaluate(-1, pt2, thetad2, phid2, trigeffnL, trigeffnH, trigeffnA);
142             thetas2 *= TMath::Pi()/180.;
143             phis2 *= TMath::Pi()/180.;
144             if (TMath::Tan(thetas2/2.) > 0) {
145                 etas  = -TMath::Log(TMath::Tan(thetas2/2.));
146                 etaH->Fill(etas, eta2);
147                 etadH->Fill(etas-eta2);
148             }
149             
150             dphi   = phis1 - phis2;
151             etas1  = - TMath::Log(TMath::Abs(TMath::Tan(thetas1/2.)+1e-4));
152             etas2  = - TMath::Log(TMath::Abs(TMath::Tan(thetas2/2.)+1e-4));         
153             deta   = etas1 - etas2;
154             pts1   = ps1 * TMath::Sin(thetas1);
155             pts2   = ps2 * TMath::Sin(thetas2);     
156             
157             m      =  TMath::Sqrt(2. * pts1 * pts2 * (TMath::CosH(deta) - TMath::Cos(dphi)));
158             Float_t wgt = effn * effp * accn * accp * trigeffpH * trigeffnH;
159             printf("Mass %f\n", m);
160             if (pts1 > 0. && pts2 > 0.)
161             massHS->Fill(m, wgt);
162         }
163         
164 //  Finish event
165 //      I/O
166 //      
167         stack->FinishEvent();
168         header->SetStack(stack);
169         rl->TreeE()->Fill();
170         header->SetNprimary(stack->GetNprimary());
171         header->SetNtrack(stack->GetNtrack());  
172         rl->WriteKinematics("OVERWRITE");
173     } // event loop
174     TCanvas *c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
175     massHU->Draw();
176     massHS->Draw("same");
177
178     TCanvas *c2 = new TCanvas("c2","Canvas 2",400,10,600,700);
179     etaH->Draw();
180     
181 //
182 //                         Termination
183 //  Generator
184     gener->FinishRun();
185 //  Stack
186     stack->FinishRun();
187 //  Write file
188     rl->WriteHeader("OVERWRITE");
189     gener->Write();
190     rl->Write();
191 }
192
193
194 AliGenerator*  CreateGenerator()
195 {
196     AliGenParam *gener 
197         = new AliGenParam(2000, AliGenMUONlib::kUpsilon, "");
198         
199     gener->SetMomentumRange(0,999);
200     gener->SetPtRange(0,100.);
201     gener->SetPhiRange(0., 360.);
202     gener->SetYRange(2.5,4);
203     gener->SetCutOnChild(1);
204     gener->SetChildThetaRange(2.0,9.);
205     gener->SetForceDecay(kDiMuon);
206
207     AliDecayer *decayer = new AliDecayerPythia();
208     decayer->SetForceDecay(kAll);
209     decayer->Init();
210     gener->SetDecayer(decayer);
211     gener->Init();
212     
213     gener->Draw();
214     
215     return gener;
216 }
217
218
219
220
221
222
223
224
225
226
227
228