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