First commit.
authormorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 12 Jun 2003 07:48:53 +0000 (07:48 +0000)
committermorsch <morsch@f7af4fe6-9843-0410-8265-dc069ae4e863>
Thu, 12 Jun 2003 07:48:53 +0000 (07:48 +0000)
FASTSIM/fastGenQuarkonia.C [new file with mode: 0644]

diff --git a/FASTSIM/fastGenQuarkonia.C b/FASTSIM/fastGenQuarkonia.C
new file mode 100644 (file)
index 0000000..52e7af3
--- /dev/null
@@ -0,0 +1,212 @@
+AliGenerator*  CreateGenerator();
+
+void fastGenQuarkonia(Int_t nev = 1, char* filename = "galice.root")
+{
+//  Update data base
+    AliPDG::AddParticlesToPdgDataBase();
+    
+// Create the fast tracker
+    AliFastMuonTriggerEff *trigeff = new AliFastMuonTriggerEff();
+    AliFastMuonTrackingAcc *acc = new AliFastMuonTrackingAcc();
+    AliFastMuonTrackingEff *eff = new AliFastMuonTrackingEff();
+    AliFastMuonTrackingRes *res = new AliFastMuonTrackingRes();
+    acc->SetBackground(0);
+    eff->SetBackground(0);
+    res->SetBackground(0);  
+    acc->Init(); 
+    eff->Init(); 
+    res->Init(); 
+    AliFastDetector* tracker = new AliFastDetector("Tracker", "Muon Tracker");
+    tracker->AddResponse(trigeff);
+    tracker->AddResponse(acc);
+    tracker->AddResponse(eff);
+    tracker->AddResponse(res);
+    tracker->Init();
+    tracker->Dump();
+//  Histos
+    TH1F* massHU = new TH1F("massHU", "Mass Spectrum:                ", 500, 5, 15.);
+    TH1F* massHS = new TH1F("massHS", "Mass Spectrum Smeared:        ", 500, 5, 15.);
+//
+//                        Construction
+//
+//  Output file
+    TFile*  file         = new TFile(filename, "recreate");
+//  Create stack
+    AliStack* stack      = new AliStack(10000);
+    stack->MakeTree(0, filename);
+
+//  Create Header
+    AliHeader* header    = new AliHeader();
+//  Create Header Tree
+    TTree* treeE         = new TTree("TE","Headers");
+    treeE->Branch("Header", "AliHeader", &header, 4000, 0);
+    treeE->Write();
+//
+//  Create and Initialize Generator
+    AliGenerator *gener = CreateGenerator();
+    AliPythia* pyth = AliPythia::Instance();
+    gener->SetStack(stack);
+    
+//
+//                        Event Loop
+//
+    Int_t iev;
+     
+    for (iev = 0; iev < nev; iev++) {
+
+       printf("\n \n Event number %d \n \n", iev);
+       
+//  Initialize event
+       header->Reset(0,iev);
+       stack->Reset();
+       stack->BeginEvent(iev);
+
+//  Generate event
+       gener->Generate();
+               
+//  Analysis
+       Int_t npart = stack->GetNprimary();
+//     printf("Analyse %d Particles\n", npart);
+       for (Int_t part=0; part<npart; part++) {
+           TParticle *MPart = stack->Particle(part);
+           Int_t mpart  = MPart->GetPdgCode();
+           if (mpart != 553 && mpart != 100553 && mpart != 200553) continue;
+           Int_t ch1 = MPart->GetFirstDaughter();
+           Int_t ch2 = MPart->GetLastDaughter();
+           
+           TParticle *muon1 = stack->Particle(ch1);
+           TParticle *muon2 = stack->Particle(ch2);
+           
+           Float_t theta1 = muon1->Theta();
+           Float_t thetad1 = theta1 * 180./TMath::Pi();
+           Float_t eta1   = muon1->Eta();
+           Float_t pt1    = muon1->Pt();
+           Float_t phi1   = muon1->Phi();
+           Float_t phid1  = phi1 * 180./TMath::Pi() - 180.;
+           Float_t p1     = muon1->P();
+
+           if (thetad1 > 9. || thetad1 < 2.) continue;
+
+           Float_t theta2 = muon2->Theta();
+           Float_t thetad2 = theta2 * 180./TMath::Pi();
+           Float_t eta2   = muon2->Eta();
+           Float_t pt2    = muon2->Pt();
+           Float_t phi2   = muon2->Phi();
+           Float_t phid2  = phi2 * 180./TMath::Pi() - 180.;
+           Float_t p2     = muon2->P();
+
+
+           if (thetad2 > 9. || thetad2 < 2.) continue;     
+           
+           Float_t dphi   = phi1 - phi2;
+           Float_t deta   = eta1 - eta2;
+           Float_t m      =  TMath::Sqrt(2. * pt1 * pt2 * (TMath::CosH(deta) - TMath::Cos(dphi)));
+
+           printf("Mass %f %f %f\n", m, phid1, phid2);
+           massHU->Fill(m);
+// Smear
+           // the mu+
+           Float_t thetas1, phis1, ps1, thetas2, phis2, ps2, pts1, pts2, etas1, etas2;
+           Float_t trigeffpL, trigeffpH, trigeffnL, trigeffnH;
+           
+           res->SetCharge(1);
+           eff->SetCharge(1);
+           acc->SetCharge(1);
+           res->Evaluate(p1, thetad1, phid1, ps1, thetas1, phis1);
+           Float_t effp     = eff->Evaluate(pt1, thetad1, phid1);
+           Float_t accp     = acc->Evaluate(pt1, thetad1, phid1);
+           trigeff->Evaluate(1, pt1, thetad1, phid1, trigeffpL, trigeffpH);
+           thetas1 *= TMath::Pi()/180.;
+           phis1 *= TMath::Pi()/180.;
+
+           // the mu- 
+           res->SetCharge(-1);
+           acc->SetCharge(-1);
+           eff->SetCharge(-1);
+           res->Evaluate(p2, thetad2, phid2, ps2, thetas2, phis2);
+           Float_t effn     = eff->Evaluate(pt2, thetad2, phid2);
+           Float_t accn     = acc->Evaluate(pt2, thetad2, phid2);
+           trigeff->Evaluate(-1, pt2, thetad2, phid2, trigeffnL, trigeffnH);
+           thetas2 *= TMath::Pi()/180.;
+           phis2 *= TMath::Pi()/180.;
+           dphi   = phis1 - phis2;
+           etas1  = - TMath::Log(TMath::Tan(thetas1/2.));
+           etas2  = - TMath::Log(TMath::Tan(thetas2/2.));          
+           deta   = etas1 - etas2;
+           pts1   = ps1 * TMath::Sin(thetas1);
+           pts2   = ps2 * TMath::Sin(thetas2);     
+           
+           m      =  TMath::Sqrt(2. * pts1 * pts2 * (TMath::CosH(deta) - TMath::Cos(dphi)));
+//         Float_t wgt = effn * effp * accn * accp * trigeffpH * trigeffnH;
+           Float_t wgt = accn * accp * effp * effn ;
+           printf("Mass %f\n", m);
+           if (pts1 > 0. && pts2 > 0.)
+           massHS->Fill(m, wgt);
+       }
+       
+//  Finish event
+       header->SetNprimary(stack->GetNprimary());
+       header->SetNtrack(stack->GetNtrack());  
+//      I/O
+//     
+//     stack->FinishEvent();
+//     header->SetStack(stack);
+//     treeE->Fill();
+//     (stack->TreeK())->Write(0,TObject::kOverwrite);
+    } // event loop
+    TCanvas *c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
+    massHU->Draw();
+    massHS->Draw("same");
+    
+//
+//                         Termination
+//  Generator
+    gener->FinishRun();
+//  Header
+    treeE->Write(0,TObject::kOverwrite);
+    delete treeE;   treeE = 0;
+//  Stack
+    stack->FinishRun();
+//  Write file
+    gener->Write();
+    file->Write();
+}
+
+
+AliGenerator*  CreateGenerator()
+{
+    AliGenParam *gener 
+       = new AliGenParam(200000, AliGenMUONlib::kUpsilon, "");
+        
+    gener->SetMomentumRange(0,999);
+    gener->SetPtRange(0,100.);
+    gener->SetPhiRange(-180, 180);
+//    gener->SetYRange(2.5,4);
+    gener->SetYRange(-6. , 6.);
+    gener->SetCutOnChild(0);
+    gener->SetChildThetaRange(2.0,9.);
+    gener->SetForceDecay(kDiMuon);
+
+    AliDecayer *decayer = new AliDecayerPythia();
+    decayer->SetForceDecay(kAll);
+    decayer->Init();
+    gener->SetDecayer(decayer);
+    gener->Init();
+    
+    gener->Draw();
+    
+    return gener;
+}
+
+
+
+
+
+
+
+
+
+
+
+