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