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
  10     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  11
  12     if (!file) {
  13         printf("\n Creating galice.root \n");
  14         file = new TFile("galice.root");
  15     } else {
  16         printf("\n galice.root found in file list");
  17     }
  18     file->ls();
  19 //
  20     TDatabasePDG*  DataBase = new TDatabasePDG();
  21
  22 // Get AliRun object from file or create it if not on file
  23     if (!gAlice) {
  24         gAlice = (AliRun*)(file->Get("gAlice"));
  25         if (gAlice) printf("AliRun object found on file\n");
  26         if (!gAlice) {
  27             printf("\n create new gAlice object");
  28             gAlice = new AliRun("gAlice","Alice test program");
  29         }
  30     }
  31
  32 //  Create some histograms
  33
  34     TH1F *thetaH   =  new TH1F("thetaH","Theta distribution",180,0,180);
  35     TH1F *phiH     =  new TH1F("phiH","Phi distribution"  ,180,-180,180);
  36     TH1F *etaH     =  new TH1F("etaH","Pseudorapidity",120,-12,12);
  37     TH1F *yH       =  new TH1F("yH","Rapidity distribution",120,-12,12);
  38     TH1F *eH       =  new TH1F("eH","Energy distribution",100,0,100);
  39     TH1F *eetaH    =  new TH1F("eetaH","Pseudorapidity",120,-12,12);
  40     TH1F *ptH      =  new TH1F("ptH","Pt distribution",150,0,15);
  41 //
  42 //   Loop over events
  43 //
  44
  45     for (Int_t nev=0; nev<= evNumber2; nev++) {
  46         Int_t nparticles = gAlice->GetEvent(nev);
  47         if (nev < evNumber1) continue;
  48         if (nparticles <= 0) return;
  49         TClonesArray *fPartArray = gAlice->Particles();
  50         Int_t npart = fPartArray->GetEntriesFast();
  51 //
  52 // Loop over primary particles (jpsi. upsilon, ...)
  53 //
  54         for (Int_t part=0; part<npart; part++) {
  55             TParticle *MPart = (TParticle*) fPartArray->UncheckedAt(part);
  56             Int_t mpart  = MPart->GetPdgCode();
  57             Int_t child1 = MPart->GetFirstDaughter();
  58             Int_t child2 = MPart->GetLastDaughter();
  59             Int_t mother = MPart->GetFirstMother();
  60
  61             Float_t Pt = MPart->Pt();
  62             Float_t E  = MPart->Energy();
  63             Float_t Pz = MPart->Pz();
  64             Float_t Py = MPart->Py();
  65             Float_t Px = MPart->Px();
  66             Float_t pT = TMath::Sqrt(Px*Px+Py*Py);
  67             Float_t theta = MPart->Theta();
  68             Float_t phi   = MPart->Phi()-TMath::Pi();
  69             Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
  70             Float_t y     = TMath::Log((E+Pz)/(E-Pz+1.e-13));
  71
  72             thetaH->Fill(theta*180./TMath::Pi(),1.);
  73             phiH->Fill(phi*180./TMath::Pi(),1.);
  74             etaH->Fill(eta,1.);
  75             eetaH->Fill(eta,E);
  76             yH->Fill(y,1.);
  77             eH->Fill(E,1.);
  78             ptH->Fill(pT,1.);
  79         } // primary loop
  80     }
  81
  82 //Create a canvas, set the view range, show histograms
  83     TCanvas *c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
  84     c1->Divide(2,2);
  85     c1->cd(1); ptH->Draw();
  86     c1->cd(2); etaH->Draw();
  87     c1->cd(3); yH->Draw();
  88     c1->cd(4
  89
  90
  91
  92 ); eH->Draw();
  93
  94     TCanvas *c2 = new TCanvas("c2","Canvas 1",400,10,600,700);
  95     c2->Divide(2,2);
  96     c2->cd(1); phiH->Draw();
  97     c2->cd(2); thetaH->Draw();
  98     c2->cd(3); eetaH->Draw();
  99 }
 100
 101
 102
 103
 104
 105
 106