MUON/MUONpadtest.C

   1 void MUONpadtest (Int_t evNumber1=0,Int_t evNumber2=1)
   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
  11   Int_t NpadX = 252;                 // number of pads on X
  12   Int_t NpadY = 374;                 // number of pads on Y
  13
  14     Int_t Pad[NpadX][NpadY];
  15     for (Int_t i=0;i<NpadX;i++) {
  16       for (Int_t j=0;j<NpadY;j++) {
  17           Pad[i][j]=0;
  18       }
  19     }
  20
  21
  22
  23 // Dynamically link some shared libs
  24
  25    if (gClassTable->GetID("AliRun") < 0) {
  26       gSystem->Load("libGeant3Dummy.so");        // a dummy version of Geant3
  27       gSystem->Load("PHOS/libPHOSdummy.so");     // the standard Alice classes
  28       gSystem->Load("libgalice.so");             // the standard Alice classes
  29    }
  30
  31 // Connect the Root Galice file containing Geometry, Kine and Hits
  32
  33    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  34    if (!file) file = new TFile("galice.root");
  35
  36 // Get AliRun object from file or create it if not on file
  37
  38    if (!gAlice) {
  39       gAlice = (AliRun*)file->Get("gAlice");
  40       if (gAlice) printf("AliRun object found on file\n");
  41       if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  42    }
  43
  44 //  Create some histograms
  45
  46    Int_t xmin=-NpadX/2;
  47    Int_t xmax= NpadX/2;
  48    Int_t ymin=-NpadY/2;
  49    Int_t ymax= NpadY/2;
  50
  51    TH2F *hc = new TH2F("hc","Chamber1 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
  52    TH2F *h = new TH2F("h","Chamber1 hit distribution",100,-100,100,100,-100,100);
  53    TH1F *charge = new TH1F("charge","Charge distribution",100,0.,1000.);
  54
  55 //   Start loop over events
  56
  57    Int_t Nh=0;
  58    Int_t Nh1=0;
  59    for (int nev=0; nev<= evNumber2; nev++) {
  60       Int_t nparticles = gAlice->GetEvent(nev);
  61       cout<<"nev  "<<nev<<endl;
  62       cout<<"nparticles  "<<nparticles<<endl;
  63       if (nev < evNumber1) continue;
  64       if (nparticles <= 0) return;
  65
  66 // Get pointers to MUON detector and Hits containers
  67
  68       AliMUON *MUON  = gAlice->GetDetector("MUON");
  69       TTree *TH = gAlice->TreeH();
  70       Int_t ntracks = TH->GetEntries();
  71
  72 // Start loop on tracks in the hits containers
  73
  74       //      Int_t Nh=0;
  75       Int_t Nc=0;
  76       for (Int_t track=0; track<ntracks;track++) {
  77         //          printf("ntracks %d\n",ntracks);
  78           gAlice->ResetHits();
  79           Int_t nbytes += TH->GetEvent(track);
  80           if (MUON)  {
  81               TClonesArray *Clusters = MUON->Clusters(); // Cluster branch
  82               TClonesArray *Hits = MUON->Hits();         // Hits branch
  83               //printf("%d %d \n",Clusters,Hits);
  84           }
  85           //see hits distribution
  86           Int_t nhits = Hits->GetEntriesFast();
  87           if (nhits) Nh+=nhits;
  88           //          printf("nhits %d\n",nhits);
  89           for (Int_t hit=0;hit<nhits;hit++) {
  90               mHit = (AliMUONhit*) Hits->UncheckedAt(hit);
  91               Int_t nch  = mHit->fChamber;              // chamber number
  92               Float_t x  = mHit->fX;                    // x-pos of hit
  93               Float_t y  = mHit->fY;                    // y-pos
  94               // Fill the histograms
  95               if( nch==1) {
  96                  Float_t rhit=TMath::Sqrt(x*x+y*y);
  97                  if( rhit<= 55 ) Nh1+=nhits;
  98                   h->Fill(x,y,(float) 1);
  99               }
 100           }
 101           //    see signal distribution
 102           Int_t nclust = Clusters->GetEntriesFast();
 103           //       printf("nclust %d\n",nclust);
 104           if (nclust) {
 105             Nc+=nclust;
 106             for (Int_t hit=0;hit<nclust;hit++) {
 107               padHit = (AliMUONcluster*) Clusters->UncheckedAt(hit);
 108               Int_t nchamber = padHit->fChamber;     // chamber number
 109               Int_t nhit = padHit->fHitNumber;          // hit number
 110               Int_t qtot = padHit->fQ;                // charge
 111               Int_t ipx  = padHit->fPadX;               // pad number on X
 112               Int_t ipy  = padHit->fPadY;               // pad number on Y
 113               Int_t iqpad  = padHit->fQpad;           // charge per pad
 114               Int_t rpad  = padHit->fRpad;            // R-position of pad
 115               if (qtot > 0 && nchamber==1){
 116                  charge->Fill(qtot,(float) 1);
 117               }
 118               if(rpad <= 55 && nchamber==1) {
 119               //              if(nchamber==1) {
 120                   Pad[ipx+126][ipy+187]+=(iqpad);
 121                   hc->Fill(ipx,ipy,(float) iqpad);
 122               }
 123             }
 124           }
 125        }
 126       //      cout<<"Nh  "<<Nh<<endl;
 127       cout<<"Nc  "<<Nc<<endl;
 128    }
 129 //Create a canvas, set the view range, show histograms
 130    TCanvas *c1 = new TCanvas("c1","Alice MUON pad hits",400,10,600,700);
 131    hc->SetXTitle("ix (npads)");
 132    hc->Draw();
 133    TCanvas *c2 = new TCanvas("c2","Alice MUON hits",400,10,600,700);
 134    h->SetFillColor(42);
 135    h->SetXTitle("x (cm)");
 136    h->Draw();
 137    TCanvas *c3 = new TCanvas("c3","Charge distribution",400,10,600,700);
 138    charge->SetFillColor(42);
 139    charge->SetXTitle("ADC units");
 140    charge->Draw();
 141
 142
 143    // calculate the number of pads which give a signal
 144
 145
 146     Int_t Np=0;
 147     for (Int_t i=0;i< NpadX;i++) {
 148       for (Int_t j=0;j< NpadY;j++) {
 149          if (Pad[i][j]>=6){
 150              Np+=1;
 151              //             cout<<"i, j  "<<i<<"  "<<j<<endl;
 152          }
 153       }
 154     }
 155
 156     cout<<"The total number of pads which give a signal "<<Np<<endl;
 157       cout<<"Nh  "<<Nh<<endl;
 158       cout<<"Nh1  "<<Nh1<<endl;
 159
 160 }