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="picts/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 = new TFile("galice.root","read");
  27   if (!file) {
  28     cerr<<"anal.C: No galice.root file. Please run the simulation\n";
  29     return 1;
  30   }
  31
  32   // Get AliRun object from file
  33   if (gAlice) delete gAlice;
  34   gAlice = (AliRun*)file->Get("gAlice");
  35   if (gAlice)
  36     cerr<<"anal.C: AliRun object found on file\n";
  37   else {
  38     cerr<<"anal.C: AliRun object not found on file\n";
  39     return 2;
  40   }
  41
  42   // Import the Kine and Hits Trees for the event evNumber in the file
  43   Int_t nparticles = gAlice->GetEvent(evNumber);
  44   if (nparticles <= 0) return;
  45
  46   // Get pointers to Alice detectors and Hits containers
  47   AliTPC *TPC  = gAlice->GetDetector("TPC");
  48   AliTRD *TRD  = gAlice->GetDetector("TRD");
  49   AliTRDgeometry *TRDgeometry   = TRD->GetGeometry();
  50
  51
  52   TTree *TH = gAlice->TreeH();
  53   Int_t ntracks    = TH->GetEntries();
  54
  55   // Create histograms
  56   TH1F *hSectors = new TH1F("hSectors","Number of tracks hits per sector",72,-0.5,71.5);
  57   TH1F *hTRD     = new TH1F("hTRD","Number of planes crossed per track in TRD",6,-0.5,5.5);
  58   TH1F *hTRDprim = new TH1F("hTRDprim","Number of planes crossed per primary track in TRD",6,-0.5,5.5);
  59
  60   // Start loop on tracks in the hits containers
  61   for (Int_t track=0; track<ntracks;track++) {
  62
  63     // ======>Histogram TPC
  64     // histogram the number of tracks per sector
  65     AliTPChit  *tpcHit;
  66      if (TPC) {
  67        for (tpcHit=(AliTPChit*)TPC->FirstHit(track);tpcHit;
  68             tpcHit=(AliTPChit*)TPC->NextHit()) {
  69          Int_t sector = tpcHit->fSector;
  70          hSectors->Fill(sector);
  71        }
  72      }
  73      // =======>Histogram TRD
  74      // histogram number of planes crossed per track
  75      // same for primary tracks only
  76      AliTRDhit  *trdHit;
  77      if (TRD) {
  78        for (trdHit=(AliTRDhit*)TRD->FirstHit(track);trdHit;
  79             trdHit=(AliTRDhit*)TRD->NextHit()) {
  80          Int_t ipart    = trdHit->Track();
  81          TParticle * particle = (TParticle*)gAlice->Particle(ipart);
  82          UShort_t detector = trdHit->GetDetector();
  83          Int_t plane = TRDgeometry->GetPlane(detector);
  84          hTRD->Fill(plane);
  85          if (particle->GetFirstMother() < 0) hTRDprim->Fill(plane);
  86        }
  87      }
  88    }
  89
  90   //Create a canvas, set the view range, show histograms
  91    TCanvas *c1 = new TCanvas("c1","Alice TPC and ITS hits",400,10,600,700);
  92    c1->Divide(1,2);
  93    c1->cd(1);
  94    gPad->SetFillColor(33);
  95    hSectors->SetFillColor(42);
  96    hSectors->Fit("pol1");
  97    c1->cd(2);
  98    gPad->SetFillColor(33);
  99    hTRD->SetFillColor(42);
 100    hTRD->Draw();
 101    hTRDprim->SetFillColor(46);
 102    hTRDprim->Draw("same");
 103    c1->Print("anal.ps");
 104 }
 105
 106
 107
 108
 109
 110
 111
 112
 113
 114
 115