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         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     = 0.5*TMath::Log((E+Pz)/(E-Pz+1.e-13));
  71
  72             if (child1 >= 0) continue;
  73             if (mpart == kPi0 || mpart == kGamma) continue;
  74             Float_t wgt = 1./(Float_t ((evNumber2-evNumber1)+1.));
  75             thetaH->Fill(theta*180./TMath::Pi(),wgt);
  76             phiH->Fill(phi*180./TMath::Pi(),wgt);
  77             etaH->Fill(eta,5.*wgt);
  78             eetaH->Fill(eta,E);
  79             yH->Fill(y,5.*wgt);
  80             eH->Fill(E,wgt);
  81             ptH->Fill(pT,wgt);
  82         } // primary loop
  83     }
  84
  85 //Create a canvas, set the view range, show histograms
  86     TCanvas *c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
  87     c1->Divide(2,2);
  88     c1->cd(1); ptH->Draw();
  89     c1->cd(2); etaH->Draw();
  90     c1->cd(3); yH->Draw();
  91     c1->cd(4
  92
  93
  94
  95 ); eH->Draw();
  96
  97     TCanvas *c2 = new TCanvas("c2","Canvas 1",400,10,600,700);
  98     c2->Divide(2,2);
  99     c2->cd(1); phiH->Draw();
 100     c2->cd(2); thetaH->Draw();
 101     c2->cd(3); eetaH->Draw();
 102     Float_t i0,i1,i2,i3,i4,i5;
 103     Float_t ne = Float_t ((evNumber2-evNumber1)+1.);
 104
 105     i0=eetaH->Integral(01,25)/1000./ne;
 106     i1=eetaH->Integral(25,40)/1000./ne;
 107     i2=eetaH->Integral(40,50)/1000./ne;
 108     i3=eetaH->Integral(50,60)/1000./ne;
 109     i4=eetaH->Integral(60,70)/1000./ne;
 110     i5=eetaH->Integral(70,120)/1000./ne;
 111
 112     printf("  0 < eta < 2.5 %f \n",i0);
 113     printf("2.5 < eta < 4.0 %f \n",i1);
 114     printf("4.0 < eta < 5.0 %f \n",i2);
 115     printf("5.0 < eta < 6.0 %f \n",i3);
 116     printf("6.0 < eta < 7.0 %f \n",i4);
 117     printf("7.0 < eta <12.0 %f \n",i5);
 118
 119 }
 120
 121
 122
 123
 124
 125
 126