Extra header added to the list
[u/mrichter/AliRoot.git] / FASTSIM / readPrimaries.C
CommitLineData
87d2b2fc 1void readPrimaries(Int_t evNumber1=0, Int_t evNumber2=0)
c5f5492e 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 AliRunLoader* rl = AliRunLoader::Open("galice.root");
11//
12 TDatabasePDG* DataBase = new TDatabasePDG();
13
14// Create some histograms
15
16 TH1F *thetaH = new TH1F("thetaH","Theta distribution",180,0,180);
17 TH1F *phiH = new TH1F("phiH","Phi distribution" ,180,-180,180);
18 TH1F *etaH = new TH1F("etaH","Pseudorapidity",120,-12,12);
19 TH1F *yH = new TH1F("yH","Rapidity distribution",120,-12,12);
20 TH1F *eH = new TH1F("eH","Energy distribution",100,0,100);
21 TH1F *eetaH = new TH1F("eetaH","Pseudorapidity",120,0,12);
22 TH1F *ptH = new TH1F("ptH","Pt distribution",150,0,15);
23//
24// Loop over events
25//
f8ef9dd2 26 rl->LoadKinematics();
27 rl->LoadHeader();
c5f5492e 28 for (Int_t nev=0; nev<= evNumber2; nev++) {
29 rl->GetEvent(nev);
c5f5492e 30 AliStack* stack = rl->Stack();
31 Int_t npart = stack->GetNprimary();
32 if (nev < evNumber1) continue;
33//
34// Loop over primary particles (jpsi. upsilon, ...)
35//
36
37 for (Int_t part=0; part<npart; part++) {
38 TParticle *MPart = stack->Particle(part);
39 Int_t mpart = MPart->GetPdgCode();
40 Int_t child1 = MPart->GetFirstDaughter();
41 Int_t child2 = MPart->GetLastDaughter();
42 Int_t mother = MPart->GetFirstMother();
43 printf("\n %d %d %d %d ", mpart, child1, child2, mother);
44
45 Float_t Pt = MPart->Pt();
46 Float_t E = MPart->Energy();
47 Float_t Pz = MPart->Pz();
48 Float_t Py = MPart->Py();
49 Float_t Px = MPart->Px();
50 Float_t pT = TMath::Sqrt(Px*Px+Py*Py);
51 Float_t theta = MPart->Theta();
52 Float_t phi = MPart->Phi()-TMath::Pi();
53 Float_t eta = -TMath::Log(TMath::Tan(theta/2.));
54 Float_t y = 0.5*TMath::Log((E+Pz+1.e-13)/(E-Pz+1.e-13));
55
56 if (child1 >= 0) continue;
57 if (mpart == kPi0 || mpart == kGamma) continue;
58 Float_t wgt = 1./(Float_t ((evNumber2-evNumber1)+1.));
59 thetaH->Fill(theta*180./TMath::Pi(),wgt);
60 phiH->Fill(phi*180./TMath::Pi(),wgt);
61 etaH->Fill(eta,5.*wgt);
62 eetaH->Fill(eta,E);
63 yH->Fill(y,5.*wgt);
64 eH->Fill(E,wgt);
65 ptH->Fill(pT,wgt);
66 } // primary loop
67 }
68
69//Create a canvas, set the view range, show histograms
70 TCanvas *c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
71 c1->Divide(2,2);
72 c1->cd(1); ptH->Draw();
73 c1->cd(2); etaH->Draw();
74 c1->cd(3); yH->Draw();
75 c1->cd(4
76
77
78
79); eH->Draw();
80
81 TCanvas *c2 = new TCanvas("c2","Canvas 1",400,10,600,700);
82 c2->Divide(2,2);
83 c2->cd(1); phiH->Draw();
84 c2->cd(2); thetaH->Draw();
85 c2->cd(3); eetaH->Draw();
86 Float_t i0,i1,i2,i3,i4,i5;
87 Float_t ne = Float_t ((evNumber2-evNumber1)+1.);
88
89 i0=eetaH->Integral(01,25)/1000./ne;
90 i1=eetaH->Integral(25,40)/1000./ne;
91 i2=eetaH->Integral(40,50)/1000./ne;
92 i3=eetaH->Integral(50,60)/1000./ne;
93 i4=eetaH->Integral(60,70)/1000./ne;
94 i5=eetaH->Integral(70,120)/1000./ne;
95
96 printf(" 0 < eta < 2.5 %f \n",i0);
97 printf("2.5 < eta < 4.0 %f \n",i1);
98 printf("4.0 < eta < 5.0 %f \n",i2);
99 printf("5.0 < eta < 6.0 %f \n",i3);
100 printf("6.0 < eta < 7.0 %f \n",i4);
101 printf("7.0 < eta <12.0 %f \n",i5);
102
103}
104
105
106
107
108
109
110