Works with newio now.

[u/mrichter/AliRoot.git] / EVGEN / TestPrimaries.C

void TestPrimaries(Int_t evNumber1=0, Int_t evNumber2=0) 
{
// Dynamically link some shared libs                    
    if (gClassTable->GetID("AliRun") < 0) {
	gROOT->LoadMacro("loadlibs.C");
	loadlibs();
    }
// Connect the Root Galice file containing Geometry, Kine and Hits
    AliRunLoader* rl = AliRunLoader::Open("galice.root");
//
    TDatabasePDG*  DataBase = TDatabasePDG::Instance();
    
//  Create some histograms

    TH1F *thetaH   =  new TH1F("thetaH","Theta distribution",180,0,180);
    TH1F *phiH     =  new TH1F("phiH","Phi distribution"  ,180,-180,180);
    TH1F *etaH     =  new TH1F("etaH","Pseudorapidity",120,-12,12);
    TH1F *yH       =  new TH1F("yH","Rapidity distribution",120,-12,12);
    TH1F *eH       =  new TH1F("eH","Energy distribution",100,0,100);
    TH1F *eetaH    =  new TH1F("eetaH","Pseudorapidity",120,0,12);
    TH1F *ptH      =  new TH1F("ptH","Pt distribution",150,0,15);
//
//   Loop over events 
//
    rl->LoadKinematics();
    rl->LoadHeader();    
    
    for (Int_t nev=0; nev<= evNumber2; nev++) {
	rl->GetEvent(nev);
	AliStack* stack = rl->Stack();
	Int_t npart = stack->GetNprimary();
	if (nev < evNumber1) continue;
//
// Loop over primary particles (jpsi. upsilon, ...)
//       
	
	for (Int_t part=0; part<npart; part++) {
	    TParticle *MPart = stack->Particle(part);
	    Int_t mpart  = MPart->GetPdgCode();
	    Int_t child1 = MPart->GetFirstDaughter();
	    Int_t child2 = MPart->GetLastDaughter();	
	    Int_t mother = MPart->GetFirstMother();	   
	    Float_t Pt = MPart->Pt();
	    Float_t E  = MPart->Energy();
	    Float_t Pz = MPart->Pz();
	    Float_t Py = MPart->Py();
	    Float_t Px = MPart->Px();
	    Float_t pT = TMath::Sqrt(Px*Px+Py*Py);
	    Float_t theta = MPart->Theta();
	    Float_t phi   = MPart->Phi()-TMath::Pi();
	    Float_t eta   = -TMath::Log(TMath::Tan(theta/2.));
	    Float_t y     = 0.5*TMath::Log((E+Pz+1.e-13)/(E-Pz+1.e-13));
	    
	    if (child1 >= 0) continue;
	    if (DataBase->GetParticle(mpart)->Charge() == 0) continue;
	    Float_t wgt = 1./(Float_t ((evNumber2-evNumber1)+1.));
	    thetaH->Fill(theta*180./TMath::Pi(),wgt);
	    phiH->Fill(phi*180./TMath::Pi(),wgt);
	    etaH->Fill(eta,5.*wgt);    
	    eetaH->Fill(eta,E);    
	    yH->Fill(y,5.*wgt);      
	    eH->Fill(E,wgt);      
	    ptH->Fill(pT,wgt);     
	} // primary loop
    }
    
//Create a canvas, set the view range, show histograms
    TCanvas *c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
    c1->Divide(2,2);
    c1->cd(1); ptH->Draw();
    c1->cd(2); etaH->Draw();
    c1->cd(3); yH->Draw();
    c1->cd(4); eH->Draw();

    TCanvas *c2 = new TCanvas("c2","Canvas 1",400,10,600,700);
    c2->Divide(2,2);
    c2->cd(1); phiH->Draw();
    c2->cd(2); thetaH->Draw();
    c2->cd(3); eetaH->Draw();

    TFile* f = new TFile("kineH.root", "update");
    ptH->Write();
    etaH->Write();
    yH->Write();
    phiH->Write();
    f->Close();
    

}