Updated.
[u/mrichter/AliRoot.git] / 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 }