MUON/MUONtest.C

   1 void MUONtestabso (Int_t evNumber1=0,Int_t evNumber2=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 do some analysis.
   7 //
   8 /////////////////////////////////////////////////////////////////////////
   9
  10 // Dynamically link some shared libs
  11
  12     if (gClassTable->GetID("AliRun") < 0) {
  13         gSystem->Load("$ALITOP/cern.so/lib/libpdfDUMMY.so");
  14         gSystem->Load("$ALITOP/cern.so/lib/libPythia.so");
  15         gSystem->Load("$ROOTSYS/lib/libEG.so");
  16         gSystem->Load("$ROOTSYS/lib/libEGPythia.so");
  17         gSystem->Load("libGeant3Dummy.so");    //a dummy version of Geant3
  18         gSystem->Load("PHOS/libPHOSdummy.so"); //the standard Alice classes
  19         gSystem->Load("libgalice.so");         // the standard Alice classes
  20     }
  21 // Connect the Root Galice file containing Geometry, Kine and Hits
  22
  23     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  24
  25
  26     if (!file) {
  27         printf("\n Creating galice.root \n");
  28         file = new TFile("galice.root");
  29     } else {
  30         printf("\n galice.root found in file list");
  31     }
  32     file->ls();
  33
  34 // Get AliRun object from file or create it if not on file
  35     if (!gAlice) {
  36         gAlice = (AliRun*)(file->Get("gAlice"));
  37         if (gAlice) printf("AliRun object found on file\n");
  38         if (!gAlice) {
  39             printf("\n create new gAlice object");
  40             gAlice = new AliRun("gAlice","Alice test program");
  41         }
  42     }
  43
  44 //  Create some histograms
  45
  46
  47    TH1F *zv = new TH1F("zv","z vertex"
  48                            ,100,468.,503.);
  49    TH1F *xv = new TH1F("xv","x vertex"
  50                            ,100,0.,50.);
  51    TH1F *yv  = new TH1F("yv","y vertex"
  52                            ,100,0.,50.);
  53
  54    TH1F *ip  = new TH1F("ip","ipart"
  55                            ,100,0.,50.);
  56
  57    AliMUONchamber*  iChamber;
  58 //
  59 //   Loop over events
  60 //
  61    Int_t Nh=0;
  62    Int_t Nh1=0;
  63    for (Int_t nev=0; nev<= evNumber2; nev++) {
  64        cout << "nev         " << nev <<endl;
  65        Int_t nparticles = gAlice->GetEvent(nev);
  66        cout << "nparticles  " << nparticles <<endl;
  67        if (nev < evNumber1) continue;
  68        if (nparticles <= 0) return;
  69
  70        AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
  71            printf("\n nev %d \n ", nev);
  72
  73        TTree *TH = gAlice->TreeH();
  74        Int_t ntracks = TH->GetEntries();
  75 //
  76 //   Loop over tracks
  77 //
  78        for (Int_t track=0; track<ntracks;track++) {
  79
  80            gAlice->ResetHits();
  81            Int_t nbytes += TH->GetEvent(track);
  82            if (MUON)  {
  83                for(AliMUONhit* mHit=(AliMUONhit*)MUON->FirstHit(-1);
  84                    mHit;
  85                    mHit=(AliMUONhit*)MUON->NextHit())
  86                {
  87                    Int_t   nch   = mHit->fChamber;  // chamber number
  88                    Float_t x     = mHit->fX;        // x-pos of hit
  89                    Float_t y     = mHit->fY;        // y-pos
  90                    Float_t z     = mHit->fZ;        // y-pos
  91
  92                 if (nch > 1) continue;
  93
  94                 Int_t ipart = mHit->fParticle;
  95                 TClonesArray *fPartArray = gAlice->Particles();
  96                 TParticle *Part;
  97                 Int_t ftrack = mHit->fTrack;
  98                 Int_t id = ((TParticle*) fPartArray->UncheckedAt(ftrack))->GetPdgCode();
  99                   ip->Fill((float)ipart,(float) 1);
 100
 101                 if (ipart != 3) continue;
 102                   Float_t xv = Part->Vx();      //  vertex
 103                   Float_t yv = Part->Vy();      //  vertex
 104                   Float_t zv  = Part->Vz();      // z vertex
 105                   xv->Fill(xv,(float) 1);
 106                   yv->Fill(yv,(float) 1);
 107                   zv->Fill(zv,(float) 1);
 108                }
 109
 110            }
 111        }
 112    }
 113
 114
 115 //Create a canvas, set the view range, show histograms
 116    TCanvas *c1 = new TCanvas("c1","Charge and Residuals",400,10,600,700);
 117    pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
 118    pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
 119    pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
 120    pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
 121    pad11->SetFillColor(11);
 122    pad12->SetFillColor(11);
 123    pad13->SetFillColor(11);
 124    pad14->SetFillColor(11);
 125    pad11->Draw();
 126    pad12->Draw();
 127    pad13->Draw();
 128    pad14->Draw();
 129
 130    pad11->cd();
 131    ip->SetFillColor(42);
 132    ip->SetXTitle("ipart");
 133    ip->Draw();
 134
 135    pad12->cd();
 136    xv->SetFillColor(42);
 137    xv->SetXTitle("xvert");
 138    xv->Draw();
 139
 140    pad13->cd();
 141    yv->SetFillColor(42);
 142    yv->SetXTitle("yvert");
 143    yv->Draw();
 144
 145    pad14->cd();
 146    zv->SetFillColor(42);
 147    zv->SetXTitle("zvert");
 148    zv->Draw();
 149
 150 }