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   const Int_t NpadX = 252;                 // number of pads on X
  11   const Int_t NpadY = 374;                 // number of pads on Y
  12
  13 // Dynamically link some shared libs
  14
  15    if (gClassTable->GetID("AliRun") < 0) {
  16       gROOT->LoadMacro("loadlibs.C");
  17       loadlibs();
  18    }
  19
  20    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  21    if (!file) file = new TFile("galice.root");
  22
  23 // Get AliRun object from file or create it if not on file
  24
  25    if (!gAlice) {
  26       gAlice = (AliRun*)file->Get("gAlice");
  27       if (gAlice) printf("AliRun object found on file\n");
  28       if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  29    }
  30
  31 //  Create some histograms
  32
  33    Int_t xmin=-NpadX/2;
  34    Int_t xmax= NpadX/2;
  35    Int_t ymin=-NpadY/2;
  36    Int_t ymax= NpadY/2;
  37
  38    TH2F *hc = new TH2F("hc","Chamber1 signal distribution",NpadX,xmin,xmax,NpadY,ymin,ymax);
  39    TH2F *h = new TH2F("h","Chamber1 hit distribution",100,-100,100,100,-100,100);
  40    TH1F *charge = new TH1F("charge","Charge distribution",100,0.,1000.);
  41
  42 //   Start loop over events
  43
  44    Int_t Nh=0;
  45    Int_t Nh1=0;
  46    AliMUON *pMUON  = (AliDetector*) gAlice->GetDetector("MUON");
  47    for (int nev=0; nev<= evNumber2; nev++) {
  48       Int_t nparticles = gAlice->GetEvent(nev);
  49       cout<<"nev  "<<nev<<endl;
  50       cout<<"nparticles  "<<nparticles<<endl;
  51       if (nev < evNumber1) continue;
  52       if (nparticles <= 0) return;
  53
  54 // Get pointers to MUON detector and Hits containers
  55       TTree *TH = gAlice->TreeH();
  56       Int_t ntracks = TH->GetEntries();
  57
  58 // Start loop on tracks in the hits containers
  59
  60       //      Int_t Nh=0;
  61       Int_t Nc=0;
  62       for (Int_t track=0; track<ntracks;track++) {
  63 //        printf("ntracks %d\n",ntracks);
  64           gAlice->ResetHits();
  65           Int_t nbytes += TH->GetEvent(track);
  66           if (pMUON)  {
  67               TClonesArray *Clusters = pMUON->PadHits();  // Cluster branch
  68               TClonesArray *Hits = pMUON->Hits();         // Hits branch
  69 //            printf("%d %d \n",Clusters,Hits);
  70           }
  71           //see hits distribution
  72           Int_t nhits = Hits->GetEntriesFast();
  73           if (nhits) Nh+=nhits;
  74 //        printf("nhits %d\n",nhits);
  75           for (Int_t hit=0;hit<nhits;hit++) {
  76 //            printf("hit# %d\n",hit);
  77               mHit = (AliMUONHit*) Hits->UncheckedAt(hit);
  78               Int_t nch  = mHit->Chamber();              // chamber number
  79               Float_t x  = mHit->X();                    // x-pos of hit
  80               Float_t y  = mHit->Y();                    // y-pos
  81               // Fill the histograms
  82               if( nch==1) {
  83                   Float_t rhit=TMath::Sqrt(x*x+y*y);
  84                   if( rhit<= 55 ) Nh1+=nhits;
  85                   h->Fill(x,y,(float) 1);
  86               }
  87
  88
  89           //    see signal distribution
  90               for (AliMUONPadHit* mPad =
  91                        (AliMUONPadHit*)pMUON->FirstPad(mHit,Clusters);
  92                    mPad;
  93                    mPad = (AliMUONPadHit*)pMUON->NextPad(Clusters))
  94               {
  95                   Int_t nhit   = mPad->HitNumber();          // hit number
  96                   Int_t qtot   = mPad->Q();                  // charge
  97                   Int_t ipx    = mPad->PadX();               // pad number on X
  98                   Int_t ipy    = mPad->PadY();               // pad number on Y
  99                   Int_t iqpad  = mPad->QPad();               // charge per pad
 100 //                printf("%d %d %d\n",ipx,ipy,iqpad);
 101                   if (qtot > 0 && nch == 1){
 102                       charge->Fill(float(iqpad),(float) 1);
 103                   }
 104                   if(nch == 1) {
 105                       hc->Fill(ipx,ipy,(float) iqpad);
 106                   } // if rpad
 107               } // padhits
 108           } // hitloops
 109       }
 110       //      cout<<"Nh  "<<Nh<<endl;
 111    }
 112 //Create a canvas, set the view range, show histograms
 113    TCanvas *c1 = new TCanvas("c1","Alice MUON pad hits",400,10,600,700);
 114    hc->SetXTitle("ix (npads)");
 115    hc->Draw();
 116    TCanvas *c2 = new TCanvas("c2","Alice MUON hits",400,10,600,700);
 117    h->SetFillColor(42);
 118    h->SetXTitle("x (cm)");
 119    h->Draw();
 120    TCanvas *c3 = new TCanvas("c3","Charge distribution",400,10,600,700);
 121    charge->SetFillColor(42);
 122    charge->SetXTitle("ADC units");
 123    charge->Draw();
 124
 125 }