First commit.
[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 //
30 //                        Construction
31 //
32 //  Output file
33     TFile*  file         = new TFile(filename, "recreate");
34 //  Create stack
35     AliStack* stack      = new AliStack(10000);
36     stack->MakeTree(0, filename);
37
38 //  Create Header
39     AliHeader* header    = new AliHeader();
40 //  Create Header Tree
41     TTree* treeE         = new TTree("TE","Headers");
42     treeE->Branch("Header", "AliHeader", &header, 4000, 0);
43     treeE->Write();
44 //
45 //  Create and Initialize Generator
46     AliGenerator *gener = CreateGenerator();
47     AliPythia* pyth = AliPythia::Instance();
48     gener->SetStack(stack);
49     
50 //
51 //                        Event Loop
52 //
53     Int_t iev;
54      
55     for (iev = 0; iev < nev; iev++) {
56
57         printf("\n \n Event number %d \n \n", iev);
58         
59 //  Initialize event
60         header->Reset(0,iev);
61         stack->Reset();
62         stack->BeginEvent(iev);
63
64 //  Generate event
65         gener->Generate();
66                 
67 //  Analysis
68         Int_t npart = stack->GetNprimary();
69 //      printf("Analyse %d Particles\n", npart);
70         for (Int_t part=0; part<npart; part++) {
71             TParticle *MPart = stack->Particle(part);
72             Int_t mpart  = MPart->GetPdgCode();
73             if (mpart != 553 && mpart != 100553 && mpart != 200553) continue;
74             Int_t ch1 = MPart->GetFirstDaughter();
75             Int_t ch2 = MPart->GetLastDaughter();
76             
77             TParticle *muon1 = stack->Particle(ch1);
78             TParticle *muon2 = stack->Particle(ch2);
79             
80             Float_t theta1 = muon1->Theta();
81             Float_t thetad1 = theta1 * 180./TMath::Pi();
82             Float_t eta1   = muon1->Eta();
83             Float_t pt1    = muon1->Pt();
84             Float_t phi1   = muon1->Phi();
85             Float_t phid1  = phi1 * 180./TMath::Pi() - 180.;
86             Float_t p1     = muon1->P();
87
88             if (thetad1 > 9. || thetad1 < 2.) continue;
89
90             Float_t theta2 = muon2->Theta();
91             Float_t thetad2 = theta2 * 180./TMath::Pi();
92             Float_t eta2   = muon2->Eta();
93             Float_t pt2    = muon2->Pt();
94             Float_t phi2   = muon2->Phi();
95             Float_t phid2  = phi2 * 180./TMath::Pi() - 180.;
96             Float_t p2     = muon2->P();
97
98
99             if (thetad2 > 9. || thetad2 < 2.) continue;     
100             
101             Float_t dphi   = phi1 - phi2;
102             Float_t deta   = eta1 - eta2;
103             Float_t m      =  TMath::Sqrt(2. * pt1 * pt2 * (TMath::CosH(deta) - TMath::Cos(dphi)));
104
105             printf("Mass %f %f %f\n", m, phid1, phid2);
106             massHU->Fill(m);
107 // Smear
108             // the mu+
109             Float_t thetas1, phis1, ps1, thetas2, phis2, ps2, pts1, pts2, etas1, etas2;
110             Float_t trigeffpL, trigeffpH, trigeffnL, trigeffnH;
111             
112             res->SetCharge(1);
113             eff->SetCharge(1);
114             acc->SetCharge(1);
115             res->Evaluate(p1, thetad1, phid1, ps1, thetas1, phis1);
116             Float_t effp     = eff->Evaluate(pt1, thetad1, phid1);
117             Float_t accp     = acc->Evaluate(pt1, thetad1, phid1);
118             trigeff->Evaluate(1, pt1, thetad1, phid1, trigeffpL, trigeffpH);
119             thetas1 *= TMath::Pi()/180.;
120             phis1 *= TMath::Pi()/180.;
121
122             // the mu- 
123             res->SetCharge(-1);
124             acc->SetCharge(-1);
125             eff->SetCharge(-1);
126             res->Evaluate(p2, thetad2, phid2, ps2, thetas2, phis2);
127             Float_t effn     = eff->Evaluate(pt2, thetad2, phid2);
128             Float_t accn     = acc->Evaluate(pt2, thetad2, phid2);
129             trigeff->Evaluate(-1, pt2, thetad2, phid2, trigeffnL, trigeffnH);
130             thetas2 *= TMath::Pi()/180.;
131             phis2 *= TMath::Pi()/180.;
132  
133             dphi   = phis1 - phis2;
134             etas1  = - TMath::Log(TMath::Tan(thetas1/2.));
135             etas2  = - TMath::Log(TMath::Tan(thetas2/2.));          
136             deta   = etas1 - etas2;
137             pts1   = ps1 * TMath::Sin(thetas1);
138             pts2   = ps2 * TMath::Sin(thetas2);     
139             
140             m      =  TMath::Sqrt(2. * pts1 * pts2 * (TMath::CosH(deta) - TMath::Cos(dphi)));
141 //          Float_t wgt = effn * effp * accn * accp * trigeffpH * trigeffnH;
142             Float_t wgt = accn * accp * effp * effn ;
143             printf("Mass %f\n", m);
144             if (pts1 > 0. && pts2 > 0.)
145             massHS->Fill(m, wgt);
146         }
147         
148 //  Finish event
149         header->SetNprimary(stack->GetNprimary());
150         header->SetNtrack(stack->GetNtrack());  
151 //      I/O
152 //      
153 //      stack->FinishEvent();
154 //      header->SetStack(stack);
155 //      treeE->Fill();
156 //      (stack->TreeK())->Write(0,TObject::kOverwrite);
157     } // event loop
158     TCanvas *c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
159     massHU->Draw();
160     massHS->Draw("same");
161     
162 //
163 //                         Termination
164 //  Generator
165     gener->FinishRun();
166 //  Header
167     treeE->Write(0,TObject::kOverwrite);
168     delete treeE;   treeE = 0;
169 //  Stack
170     stack->FinishRun();
171 //  Write file
172     gener->Write();
173     file->Write();
174 }
175
176
177 AliGenerator*  CreateGenerator()
178 {
179     AliGenParam *gener 
180         = new AliGenParam(200000, AliGenMUONlib::kUpsilon, "");
181         
182     gener->SetMomentumRange(0,999);
183     gener->SetPtRange(0,100.);
184     gener->SetPhiRange(-180, 180);
185 //    gener->SetYRange(2.5,4);
186     gener->SetYRange(-6. , 6.);
187     gener->SetCutOnChild(0);
188     gener->SetChildThetaRange(2.0,9.);
189     gener->SetForceDecay(kDiMuon);
190
191     AliDecayer *decayer = new AliDecayerPythia();
192     decayer->SetForceDecay(kAll);
193     decayer->Init();
194     gener->SetDecayer(decayer);
195     gener->Init();
196     
197     gener->Draw();
198     
199     return gener;
200 }
201
202
203
204
205
206
207
208
209
210
211
212