EVGEN/TestPrimaries.C

   1 void TestPrimaries(Int_t evNumber1=0, Int_t evNumber2=0)
   2 {
   3 // Dynamically link some shared libs
   4     if (gClassTable->GetID("AliRun") < 0) {
   5         gROOT->LoadMacro("loadlibs.C");
   6         loadlibs();
   7     }
   8 // Connect the Root Galice file containing Geometry, Kine and Hits
   9     AliRunLoader* rl = AliRunLoader::Open("galice.root");
  10 //
  11     TDatabasePDG*  DataBase = TDatabasePDG::Instance();
  12
  13 //  Create some histograms
  14
  15     TH1F *thetaH   =  new TH1F("thetaH","Theta distribution",180,0,180);
  16     thetaH->Sumw2();
  17     TH1F *phiH     =  new TH1F("phiH","Phi distribution"  ,180,-180,180);
  18     phiH->Sumw2();
  19     TH1F *phiH1    =  new TH1F("phiH1","Phi distribution"  ,180,-180,180);
  20     phiH1->Sumw2();
  21     TH1F *etaH     =  new TH1F("etaH","Pseudorapidity",120,-12,12);
  22     etaH->Sumw2();
  23     TH1F *yH       =  new TH1F("yH","Rapidity distribution",120,-12,12);
  24     yH->Sumw2();
  25     TH1F *eH       =  new TH1F("eH","Energy distribution",100,0,100);
  26     eH->Sumw2();
  27     TH1F *eetaH    =  new TH1F("eetaH","Pseudorapidity",120,0,12);
  28     eetaH->Sumw2();
  29     TH1F *ptH      =  new TH1F("ptH","Pt distribution",150,0,15);
  30     ptH->Sumw2();
  31 //
  32 //   Loop over events
  33 //
  34     rl->LoadKinematics();
  35     rl->LoadHeader();
  36
  37     for (Int_t nev=0; nev<= evNumber2; nev++) {
  38         rl->GetEvent(nev);
  39         AliStack* stack = rl->Stack();
  40         Int_t npart = stack->GetNprimary();
  41         if (nev < evNumber1) continue;
  42 //
  43 // Loop over primary particles (jpsi. upsilon, ...)
  44 //
  45
  46         for (Int_t part=0; part<npart; part++) {
  47             TParticle *MPart = stack->Particle(part);
  48             Int_t mpart  = MPart->GetPdgCode();
  49             Int_t child1 = MPart->GetFirstDaughter();
  50             Int_t child2 = MPart->GetLastDaughter();
  51             Int_t mother = MPart->GetFirstMother();
  52             Float_t Pt = MPart->Pt();
  53             Float_t E  = MPart->Energy();
  54             Float_t Pz = MPart->Pz();
  55             Float_t Py = MPart->Py();
  56             Float_t Px = MPart->Px();
  57             Float_t pT = TMath::Sqrt(Px*Px+Py*Py);
  58             Float_t theta = MPart->Theta();
  59             Float_t phi   = MPart->Phi()-TMath::Pi();
  60             Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
  61             Float_t y     = 0.5*TMath::Log((E+Pz+1.e-13)/(E-Pz+1.e-13));
  62
  63             if (child1 >= 0) continue;
  64             if (DataBase->GetParticle(mpart)->Charge() == 0) continue;
  65             Float_t wgt = 1./(Float_t ((evNumber2-evNumber1)+1.));
  66             thetaH->Fill(theta*180./TMath::Pi(),wgt);
  67             if (pT<2)
  68               phiH->Fill(phi*180./TMath::Pi(),wgt);
  69             else
  70               phiH1->Fill(phi*180./TMath::Pi(),wgt);
  71             etaH->Fill(eta,5.*wgt);
  72             eetaH->Fill(eta,E);
  73             yH->Fill(y,5.*wgt);
  74             eH->Fill(E,wgt);
  75             ptH->Fill(pT,wgt);
  76         } // primary loop
  77     }
  78
  79 //Create a canvas, set the view range, show histograms
  80     TCanvas *c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
  81     c1->Divide(2,2);
  82     c1->cd(1); ptH->Draw();
  83     c1->cd(2); etaH->Draw();
  84     c1->cd(3); yH->Draw();
  85     c1->cd(4); eH->Draw();
  86
  87     TCanvas *c2 = new TCanvas("c2","Canvas 1",400,10,600,700);
  88     c2->Divide(2,2);
  89     c2->cd(1); phiH->Draw();
  90     c2->cd(2); phiH1->Draw();
  91     c2->cd(3); thetaH->Draw();
  92     c2->cd(4); eetaH->Draw();
  93
  94     TFile* f = new TFile("kineH.root", "update");
  95     ptH->Write();
  96     etaH->Write();
  97     yH->Write();
  98     phiH->Write();
  99     phiH1->Write();
 100     f->Close();
 101
 102
 103 }
 104
 105
 106
 107
 108
 109
 110