Add const to SolenoidField() method for correct overwriting from the parent class
[u/mrichter/AliRoot.git] / 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+1.e-13)/(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