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