Test macro for generator output.
[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,-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