MUON/MUONtestabso.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       gROOT->LoadMacro("loadlibs.C");
  14       loadlibs();
  15    }
  16
  17 // Connect the Root Galice file containing Geometry, Kine and Hits
  18
  19     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  20
  21
  22     if (!file) {
  23         printf("\n Creating galice.root \n");
  24         file = new TFile("galice.root");
  25     } else {
  26         printf("\n galice.root found in file list");
  27     }
  28     file->ls();
  29
  30 // Get AliRun object from file or create it if not on file
  31     if (!gAlice) {
  32         gAlice = (AliRun*)(file->Get("gAlice"));
  33         if (gAlice) printf("AliRun object found on file\n");
  34         if (!gAlice) {
  35             printf("\n create new gAlice object");
  36             gAlice = new AliRun("gAlice","Alice test program");
  37         }
  38     }
  39
  40    printf ("I'm after gAlice \n");
  41
  42 //  Create some histograms
  43
  44
  45    TH1F *zv = new TH1F("zv","z vertex"
  46                            ,100,450..,506.);
  47    TH1F *xv = new TH1F("xv","x vertex"
  48                            ,140,-70.,70.);
  49    TH1F *yv  = new TH1F("yv","y vertex"
  50                            ,140,-70.,70.);
  51
  52    TH1F *ip  = new TH1F("ip","geant part"
  53                            ,50,0.,10.);
  54    TH2F *rzv  = new TH2F("rzv","R-z vert"
  55                            ,100,500.,506.,50,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                 Part = (TParticle*) fPartArray->UncheckedAt(ftrack);
  99                 //Int_t id = ((TParticle*) fPartArray->UncheckedAt(ftrack))->GetPdgCode();
 100                   ip->Fill((float)ipart);
 101
 102                 if (ipart > 3) continue;
 103
 104                   Float_t xvert = Part->Vx();      //  vertex
 105                   Float_t yvert = Part->Vy();      //  vertex
 106                   Float_t zvert  = Part->Vz();      // z vertex
 107                   xv->Fill(xvert);
 108                   yv->Fill(yvert);
 109                   zv->Fill(zvert);
 110                   Float_t rvert=TMath::Sqrt(xvert*xvert+yvert*yvert);
 111                   rzv->Fill(zvert,rvert);
 112
 113                }
 114
 115            }
 116        }
 117    }
 118
 119
 120 //Create a canvas, set the view range, show histograms
 121    TCanvas *c1 = new TCanvas("c1","Vetices from electrons and positrons",400,10,600,700);
 122    pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
 123    pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
 124    pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
 125    pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
 126    pad11->SetFillColor(11);
 127    pad12->SetFillColor(11);
 128    pad13->SetFillColor(11);
 129    pad14->SetFillColor(11);
 130    pad11->Draw();
 131    pad12->Draw();
 132    pad13->Draw();
 133    pad14->Draw();
 134
 135    pad11->cd();
 136    ip->SetFillColor(42);
 137    ip->SetXTitle("ipart");
 138    ip->Draw();
 139
 140    pad12->cd();
 141    xv->SetFillColor(42);
 142    xv->SetXTitle("xvert");
 143    xv->Draw();
 144
 145    pad13->cd();
 146    yv->SetFillColor(42);
 147    yv->SetXTitle("yvert");
 148    yv->Draw();
 149
 150    pad14->cd();
 151    zv->SetFillColor(42);
 152    zv->SetXTitle("zvert");
 153    zv->Draw();
 154
 155    TCanvas *c2 = new TCanvas("c2","R-Z vertex distribution",400,10,600,700);
 156    rzv->SetXTitle("zvert");
 157    rzv->SetYTitle("rvert");
 158    rzv->Draw();
 159 }