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,0,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
  47         Int_t nparticles = gAlice->GetEvent(nev);
  48         if (nev < evNumber1) continue;
  49         if (nparticles <= 0) return;
  50
  51         TObjArray*  parray = gAlice->Particles();
  52         Int_t npart = (gAlice->TreeK())->GetEntries();
  53 //
  54 // Loop over primary particles (jpsi. upsilon, ...)
  55 //
  56
  57         for (Int_t part=0; part<npart; part++) {
  58             TParticle *MPart = gAlice->Particle(part);
  59             Int_t mpart  = MPart->GetPdgCode();
  60             Int_t child1 = MPart->GetFirstDaughter();
  61             Int_t child2 = MPart->GetLastDaughter();
  62             Int_t mother = MPart->GetFirstMother();
  63             printf("\n  %d %d %d %d ", mpart, child1, child2, mother);
  64
  65             Float_t Pt = MPart->Pt();
  66             Float_t E  = MPart->Energy();
  67             Float_t Pz = MPart->Pz();
  68             Float_t Py = MPart->Py();
  69             Float_t Px = MPart->Px();
  70             Float_t pT = TMath::Sqrt(Px*Px+Py*Py);
  71             Float_t theta = MPart->Theta();
  72             Float_t phi   = MPart->Phi()-TMath::Pi();
  73             Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
  74             Float_t y     = 0.5*TMath::Log((E+Pz)/(E-Pz+1.e-13));
  75
  76             if (child1 >= 0) continue;
  77             if (mpart == kPi0 || mpart == kGamma) continue;
  78             Float_t wgt = 1./(Float_t ((evNumber2-evNumber1)+1.));
  79             thetaH->Fill(theta*180./TMath::Pi(),wgt);
  80             phiH->Fill(phi*180./TMath::Pi(),wgt);
  81             etaH->Fill(eta,5.*wgt);
  82             eetaH->Fill(eta,E);
  83             yH->Fill(y,5.*wgt);
  84             eH->Fill(E,wgt);
  85             ptH->Fill(pT,wgt);
  86         } // primary loop
  87     }
  88
  89 //Create a canvas, set the view range, show histograms
  90     TCanvas *c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
  91     c1->Divide(2,2);
  92     c1->cd(1); ptH->Draw();
  93     c1->cd(2); etaH->Draw();
  94     c1->cd(3); yH->Draw();
  95     c1->cd(4
  96
  97
  98
  99 ); eH->Draw();
 100
 101     TCanvas *c2 = new TCanvas("c2","Canvas 1",400,10,600,700);
 102     c2->Divide(2,2);
 103     c2->cd(1); phiH->Draw();
 104     c2->cd(2); thetaH->Draw();
 105     c2->cd(3); eetaH->Draw();
 106     Float_t i0,i1,i2,i3,i4,i5;
 107     Float_t ne = Float_t ((evNumber2-evNumber1)+1.);
 108
 109     i0=eetaH->Integral(01,25)/1000./ne;
 110     i1=eetaH->Integral(25,40)/1000./ne;
 111     i2=eetaH->Integral(40,50)/1000./ne;
 112     i3=eetaH->Integral(50,60)/1000./ne;
 113     i4=eetaH->Integral(60,70)/1000./ne;
 114     i5=eetaH->Integral(70,120)/1000./ne;
 115
 116     printf("  0 < eta < 2.5 %f \n",i0);
 117     printf("2.5 < eta < 4.0 %f \n",i1);
 118     printf("4.0 < eta < 5.0 %f \n",i2);
 119     printf("5.0 < eta < 6.0 %f \n",i3);
 120     printf("6.0 < eta < 7.0 %f \n",i4);
 121     printf("7.0 < eta <12.0 %f \n",i5);
 122
 123 }
 124
 125
 126
 127
 128
 129
 130