macros/newanal.C

   1 void newanal (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 = (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
  46 // Get pointers to Alice detectors and Hits containers
  47    AliDetector *TPC  = gAlice->GetDetector("TPC");
  48    AliDetector *TRD  = gAlice->GetDetector("TRD");
  49    TClonesArray *Particles = gAlice->Particles();
  50
  51    Int_t ntracks    = gAlice->TreeH()->GetEntries();
  52
  53    // Create histograms
  54    TH1F *hSectors = new TH1F("hSectors","Number of tracks hits per sector",75,1,76);
  55    TH1F *hTRD     = new TH1F("hTRD","Number of planes crossed per track in TRD",6,1,7);
  56    TH1F *hTRDprim = new TH1F("hTRDprim","Number of planes crossed per primary track in TRD",6,1,7);
  57
  58 // Start loop on tracks in the hits containers
  59    for (Int_t track=0; track<ntracks;track++) {
  60      if(TPC) {
  61      // ======>Histogram TPC
  62      // histogram number of tracks per sector
  63        for(AliTPChit* tpcHit=(AliTPChit*)TPC->FirstHit(track); tpcHit; tpcHit=(AliTPChit*)TPC->NextHit()) {
  64          sector = tpcHit->fPadRow;
  65          hSectors->Fill(sector);
  66        }
  67      }
  68      // =======>Histogram TRD
  69      // histogram number of planes crossed per track
  70      // same for primary tracks only
  71      if (TRD) {
  72        for(AliTRDhit* trdHit=(AliTRDhit*)TRD->FirstHit(-1); trdHit; trdHit=(AliTRDhit*)TRD->NextHit()) {
  73          ipart    = trdHit->fTrack;
  74          particle = (GParticle*)Particles->UncheckedAt(ipart);
  75          plane = trdHit->fPlane;
  76          hTRD->Fill(plane);
  77          if (particle->GetParent() < 0) hTRDprim->Fill(plane);
  78        }
  79      }
  80    }
  81
  82 //Create a canvas, set the view range, show histograms
  83    TCanvas *c1 = new TCanvas("c1","Alice TPC and ITS hits",400,10,600,700);
  84    c1->Divide(1,2);
  85    c1->cd(1);
  86    gPad->SetFillColor(33);
  87    hSectors->SetFillColor(42);
  88    hSectors->Fit("pol1");
  89    c1->cd(2);
  90    gPad->SetFillColor(33);
  91    hTRD->SetFillColor(42);
  92    hTRD->Draw();
  93    hTRDprim->SetFillColor(46);
  94    hTRDprim->Draw("same");
  95    c1->Print("anal.ps");
  96 }