macros/anal.C

   1 void anal (Int_t evNumber=0)
   2 {
   3 /////////////////////////////////////////////////////////////////////////
   4 //   This macro is a small example of a ROOT macro
   5 //   illustrating how to read the output of GALICE
   6 //   and fill some histograms.
   7 //
   8 //     Root > .L anal.C   //this loads the macro in memory
   9 //     Root > anal();     //by default process first event
  10 //     Root > anal(2);    //process third event
  11 //Begin_Html
  12 /*
  13 <img src="gif/anal.gif">
  14 */
  15 //End_Html
  16 /////////////////////////////////////////////////////////////////////////
  17
  18
  19 // Dynamically link some shared libs
  20    if (gClassTable->GetID("AliRun") < 0) {
  21       gROOT->LoadMacro("loadlibs.C");
  22       loadlibs();
  23    }
  24
  25 // Connect the Root Galice file containing Geometry, Kine and Hits
  26    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  27    if (!file) file = new TFile("galice.root");
  28
  29 // Get AliRun object from file or create it if not on file
  30    if (!gAlice) {
  31       gAlice = (AliRun*)file->Get("gAlice");
  32       if (gAlice) printf("AliRun object found on file\n");
  33       if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  34    }
  35
  36 // Import the Kine and Hits Trees for the event evNumber in the file
  37    Int_t nparticles = gAlice->GetEvent(evNumber);
  38    if (nparticles <= 0) return;
  39    Float_t x,y,z,mass,e;
  40    Int_t nbytes = 0;
  41    Int_t j,hit,ipart;
  42    Int_t nhits;
  43    Int_t sector,plane;
  44    GParticle *particle;
  45    AliTPChit  *tpcHit;
  46    AliTRDhit  *trdHit;
  47
  48 // Get pointers to Alice detectors and Hits containers
  49    AliDetector *TPC  = gAlice->GetDetector("TPC");
  50    AliDetector *TRD  = gAlice->GetDetector("TRD");
  51    TClonesArray *Particles = gAlice->Particles();
  52    if (TPC) TClonesArray *TPChits   = TPC->Hits();
  53    if (TRD) TClonesArray *TRDhits   = TRD->Hits();
  54
  55    TTree *TH = gAlice->TreeH();
  56    Int_t ntracks    = TH->GetEntries();
  57
  58    // Create histograms
  59    TH1F *hSectors = new TH1F("hSectors","Number of tracks hits per sector",75,1,76);
  60    TH1F *hTRD     = new TH1F("hTRD","Number of planes crossed per track in TRD",6,1,7);
  61    TH1F *hTRDprim = new TH1F("hTRDprim","Number of planes crossed per primary track in TRD",6,1,7);
  62
  63 // Start loop on tracks in the hits containers
  64    for (Int_t track=0; track<ntracks;track++) {
  65      gAlice->ResetHits();
  66      nbytes += TH->GetEvent(track);
  67      // ======>Histogram TPC
  68      // histogram number of tracks per sector
  69      if (TPC) {
  70        nhits = TPChits->GetEntriesFast();
  71        for (hit=0;hit<nhits;hit++) {
  72          tpcHit   = (AliTPChit*)TPChits->UncheckedAt(hit);
  73          sector = tpcHit->fPadRow;
  74          hSectors->Fill(sector);
  75        }
  76      }
  77      // =======>Histogram TRD
  78      // histogram number of planes crossed per track
  79      // same for primary tracks only
  80      if (TRD) {
  81        nhits = TRDhits->GetEntriesFast();
  82        for (hit=0;hit<nhits;hit++) {
  83          trdHit   = (AliTRDhit*)TRDhits->UncheckedAt(hit);
  84          ipart    = trdHit->fTrack;
  85          particle = (GParticle*)Particles->UncheckedAt(ipart);
  86          plane = trdHit->fPlane;
  87          hTRD->Fill(plane);
  88          if (particle->GetParent() < 0) hTRDprim->Fill(plane);
  89        }
  90      }
  91    }
  92
  93 //Create a canvas, set the view range, show histograms
  94    TCanvas *c1 = new TCanvas("c1","Alice TPC and ITS hits",400,10,600,700);
  95    c1->Divide(1,2);
  96    c1->cd(1);
  97    gPad->SetFillColor(33);
  98    hSectors->SetFillColor(42);
  99    hSectors->Fit("pol1");
 100    c1->cd(2);
 101    gPad->SetFillColor(33);
 102    hTRD->SetFillColor(42);
 103    hTRD->Draw();
 104    hTRDprim->SetFillColor(46);
 105    hTRDprim->Draw("same");
 106    c1->Print("anal.ps");
 107 }