Code from MUON-dev joined
[u/mrichter/AliRoot.git] / MUON / MUONhitrectest.C
1 #include "iostream.h"
2
3 void MUONhitrectest (Int_t evNumber1=0,Int_t evNumber2=0) 
4 {
5 /////////////////////////////////////////////////////////////////////////
6 //   This macro is a small example of a ROOT macro
7 //   illustrating how to read the output of GALICE
8 //   and do some analysis.
9 //   
10 /////////////////////////////////////////////////////////////////////////
11
12 // Dynamically link some shared libs
13
14    if (gClassTable->GetID("AliRun") < 0) {
15       gROOT->LoadMacro("loadlibs.C");
16       loadlibs();
17    }
18
19 // Connect the Root Galice file containing Geometry, Kine and Hits
20
21    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
22    if (!file) file = new TFile("galice.root");
23
24 // Get AliRun object from file or create it if not on file
25
26    if (!gAlice) {
27       gAlice = (AliRun*)file->Get("gAlice");
28       if (gAlice) printf("AliRun object found on file\n");
29       if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
30    }
31 //  Create some histograms
32
33    TH2F *h21 = new TH2F("h21","Hits",100,-100,100,100,-100,100);
34    TH2F *h22 = new TH2F("h22","CoG ",100,-100,100,100,-100,100);
35    TH1F *h1 = new TH1F("h1","Multiplicity",30,-0.5,29.5);
36    TH1F *hmult = new TH1F("hmult","Multiplicity",30,-0.5,29.5);
37    TH1F *hresx = new TH1F("hresx","Residuals",100,-1,1);
38    TH1F *hresy = new TH1F("hresy","Residuals",100,-10.,10);
39    TH1F *hresym = new TH1F("hresym","Residuals",100,-500,500);
40    TH1F *hpos = new TH1F("hnpos","Possibilities",10,-0.5,9.5);
41    TH2F *hchi1 = new TH2F("hchi1","Chi2 vs Residuals",100,0,0.2,100,-500,500);
42    TH2F *hchi2 = new TH2F("hchi2","Chi2 vs Residuals",100,0,20,100,-500,500);
43 //
44 //   Loop over events 
45 //
46    Int_t Nh=0;
47    Int_t Nh1=0;
48
49    for (int nev=0; nev<= evNumber2; nev++) {
50        Int_t nparticles = gAlice->GetEvent(nev);
51        cout << "nev         " << nev <<endl;
52        cout << "nparticles  " << nparticles <<endl;
53        if (nev < evNumber1) continue;
54        if (nparticles <= 0) return;
55        
56        TTree *TH = gAlice->TreeH();
57        Int_t ntracks = TH->GetEntries();
58        cout<<"ntracks "<<ntracks<<endl;
59
60        Int_t nbytes = 0;
61
62
63
64 // Get pointers to Alice detectors and Digits containers
65        AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
66        TClonesArray *Particles = gAlice->Particles();
67        TTree *TR = gAlice->TreeR();
68        Int_t nent=TR->GetEntries();
69        printf("Found %d entries in the tree (must be one per cathode per event! + 1empty)\n",nent);
70        if (MUON) {
71            Float_t xc[2000], yc[2000], itrc[2000];
72            TClonesArray *MUONrawclust  = MUON->RawClustAddress(0);
73            nbytes += TR->GetEvent(1);
74            Int_t nrawcl = MUONrawclust->GetEntries();
75            AliMUONRawCluster  *mRaw;
76            for (Int_t iraw=0; iraw < nrawcl; iraw++) {
77                mRaw = (AliMUONRawCluster*)MUONrawclust->UncheckedAt(iraw);
78                xc[iraw]=mRaw->fX[1];
79                yc[iraw]=mRaw->fY[0];       
80                itrc[iraw]=mRaw->fTracks[1];
81            } // cluster
82
83            gAlice->ResetHits();
84            for (Int_t itrack=0; itrack<ntracks ;itrack++) 
85            {
86                Int_t nbytes += TH->GetEvent(itrack);
87                Int_t nhit=MUON->Hits()->GetEntriesFast();
88                for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(itrack); 
89                    mHit;
90                    mHit=(AliMUONHit*)MUON->NextHit()) 
91                {
92                    Int_t   nch   = mHit->fChamber;  // chamber number
93                    Float_t x     = mHit->fX;        // x-pos of hit
94                    Float_t y     = mHit->fY;        // y-pos
95                    Int_t   ip    = mHit->fParticle;
96                    if (ip != kMuonPlus && ip != kMuonMinus) continue;
97                    if (nch != 1) continue;
98                    Float_t dmin=1.e8;
99                    Int_t imin=-1;
100 //
101 // Loop over clusters
102                    Int_t npos=0;
103                    
104                    for (Int_t i=0; i<nrawcl; i++) {
105                        Float_t dx = xc[i]-x;
106                        Float_t dy = yc[i]-y;
107                        Float_t dr  = TMath::Sqrt(dx*dx+dy*dy);
108                        if (dr<dmin) {
109                            imin=i;
110                            dmin=dr;
111                        }
112                        if (TMath::Abs(dy)< 0.05 && TMath::Abs(dx)< 0.2 ) npos++;
113                    }
114                    mRaw = (AliMUONRawCluster*)MUONrawclust->UncheckedAt(imin);
115                    Int_t   track=mRaw->fTracks[1];
116                    Float_t xrec=mRaw->fX[1];
117                    Float_t yrec=mRaw->fY[0];
118                    Float_t x_res=xrec-x;
119                    Float_t y_res=yrec-y;
120                    hresx->Fill(x_res,1.);
121                    hresy->Fill(y_res,1.);
122                    hpos->Fill(Float_t(npos), 1.);
123                    hresym->Fill(y_res*1.e4,1.);        
124                    if (npos == 0) {
125                        printf("No cluster %d %f %f %d\n", 
126                               itrack, x, y, nhit);
127                        h21->Fill(x,y,1.);
128                    }
129                    
130                    
131                    
132                }  // hits
133            } // tracks
134        }   // end if MUON
135    }   // event loop 
136    TCanvas *c1 = new TCanvas("c1","Charge and Residuals",400,10,600,700);
137    c1->Divide(2,2);
138    c1->cd(1);
139    hresx->SetFillColor(42);
140    hresx->SetXTitle("xrec-x");
141    hresx->Draw();
142
143    c1->cd(2);
144    hresy->SetFillColor(42);
145    hresy->SetXTitle("yrec-y");
146    hresy->Draw();
147
148    c1->cd(3);
149    hpos->SetFillColor(42);
150    hpos->SetXTitle("Possibilities");
151    hpos->Draw();
152
153    c1->cd(4);
154    hresym->SetFillColor(42);
155    hresym->SetXTitle("yrec-y");
156    hresym->Fit("gaus");
157    hresym->Draw();
158
159    TCanvas *c2 = new TCanvas("c2","Charge and Residuals",400,10,600,700);
160    c2->Divide(2,2);
161
162    c2->cd(1);
163    h21->SetFillColor(42);
164    h21->SetXTitle("x");
165    h21->SetYTitle("y");
166    h21->Draw();
167
168    c2->cd(2);
169    hmult->SetFillColor(42);
170    hmult->SetXTitle("multiplicity");
171    hmult->Draw();
172
173    c2->cd(3);
174    hchi1->SetFillColor(42);
175    hchi1->SetXTitle("chi2");
176    hchi1->Draw();
177
178    c2->cd(4);
179    hchi2->SetFillColor(42);
180    hchi2->SetXTitle("chi2");
181    hchi2->Draw();
182
183 }
184
185
186