Modifications needed by the HBT analysis (P.Skowronski)
[u/mrichter/AliRoot.git] / MUON / MUONhitrectest.C
1 void MUONhitrectest (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    if (!file) file = new TFile("galice.root");
21
22 // Get AliRun object from file or create it if not on file
23
24    if (!gAlice) {
25       gAlice = (AliRun*)file->Get("gAlice");
26       if (gAlice) printf("AliRun object found on file\n");
27       if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
28    }
29 //  Create some histograms
30
31    TH2F *h21 = new TH2F("h21","Hits",100,-100,100,100,-100,100);
32    TH2F *h22 = new TH2F("h22","CoG ",100,-100,100,100,-100,100);
33    TH1F *h1 = new TH1F("h1","Multiplicity",30,-0.5,29.5);
34    TH1F *hmult = new TH1F("hmult","Multiplicity",30,-0.5,29.5);
35    TH1F *hresx = new TH1F("hresx","Residuals",100,-1,1);
36    TH1F *hresy = new TH1F("hresy","Residuals",100,-10.,10);
37    TH1F *hresym = new TH1F("hresym","Residuals",100,-500,500);
38    TH1F *hpos = new TH1F("hnpos","Possibilities",10,-0.5,9.5);
39    TH2F *hchi1 = new TH2F("hchi1","Chi2 vs Residuals",100,0,0.2,100,-500,500);
40    TH2F *hchi2 = new TH2F("hchi2","Chi2 vs Residuals",100,0,20,100,-500,500);
41 //
42 //   Loop over events 
43 //
44    Int_t Nh=0;
45    Int_t Nh1=0;
46
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        TTree *TH = gAlice->TreeH();
55        Int_t ntracks = TH->GetEntries();
56        cout<<"ntracks "<<ntracks<<endl;
57
58        Int_t nbytes = 0;
59
60
61
62 // Get pointers to Alice detectors and Digits containers
63        AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
64        TClonesArray *Particles = gAlice->Particles();
65        TTree *TR = gAlice->TreeR();
66        Int_t nent=TR->GetEntries();
67
68        if (MUON) {
69            Float_t xc[2000], yc[2000], itrc[2000];
70            TClonesArray *MUONrawclust  = MUON->RawClustAddress(0);
71            nbytes += TR->GetEvent(0);
72            Int_t nrawcl = MUONrawclust->GetEntries();
73            AliMUONRawCluster  *mRaw;
74            for (Int_t iraw=0; iraw < nrawcl; iraw++) {
75                mRaw = (AliMUONRawCluster*)MUONrawclust->UncheckedAt(iraw);
76                xc[iraw]=mRaw->fX[1];
77                yc[iraw]=mRaw->fY[0];       
78                itrc[iraw]=mRaw->fTracks[1];
79            } // cluster
80
81            gAlice->ResetHits();
82            for (Int_t itrack=0; itrack<ntracks ;itrack++) 
83            {
84                Int_t nbytes += TH->GetEvent(itrack);
85                Int_t nhit=MUON->Hits()->GetEntriesFast();
86                for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(itrack); 
87                    mHit;
88                    mHit=(AliMUONHit*)MUON->NextHit()) 
89                {
90                    Int_t   nch   = mHit->Chamber();  // chamber number
91                    Float_t x     = mHit->X();        // x-pos of hit
92                    Float_t y     = mHit->Y();        // y-pos
93                    Int_t   ip    = mHit->Particle();
94                    if (ip != kMuonPlus && ip != kMuonMinus) continue;
95                    if (nch != 1) continue;
96                    Float_t dmin=1.e8;
97                    Int_t imin=-1;
98 //
99 // Loop over clusters
100                    Int_t npos=0;
101                    
102                    for (Int_t i=0; i<nrawcl; i++) {
103                        Float_t dx = xc[i]-x;
104                        Float_t dy = yc[i]-y;
105                        Float_t dr = TMath::Sqrt(dx*dx+dy*dy);
106                        if (dr<dmin) {
107                            imin=i;
108                            dmin=dr;
109                        }
110                        if (TMath::Abs(dy) < 0.05 && TMath::Abs(dx) < 0.2 ) npos++;
111                    }
112                    mRaw = (AliMUONRawCluster*)MUONrawclust->UncheckedAt(imin);
113                    Int_t   track = mRaw->fTracks[1];
114                    Float_t xrec  = mRaw->fX[1];
115                    Float_t yrec  = mRaw->fY[0];
116                    Float_t x_res = xrec-x;
117                    Float_t y_res = yrec-y;
118                    if (y_res > 0.05 ) printf("\n track %d %f\n", itrack, y_res);
119                    
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                }  // hits
130            } // tracks
131        }   // end if MUON
132    }   // event loop 
133    TCanvas *c1 = new TCanvas("c1","Charge and Residuals",400,10,600,700);
134    c1->Divide(2,2);
135    c1->cd(1);
136    hresx->SetFillColor(42);
137    hresx->SetXTitle("xrec-x");
138    hresx->Draw();
139
140    c1->cd(2);
141    hresy->SetFillColor(42);
142    hresy->SetXTitle("yrec-y");
143    hresy->Draw();
144
145    c1->cd(3);
146    hpos->SetFillColor(42);
147    hpos->SetXTitle("Possibilities");
148    hpos->Draw();
149
150    c1->cd(4);
151    hresym->SetFillColor(42);
152    hresym->SetXTitle("yrec-y");
153    hresym->Fit("gaus");
154    hresym->Draw();
155
156    TCanvas *c2 = new TCanvas("c2","Charge and Residuals",400,10,600,700);
157    c2->Divide(2,2);
158
159    c2->cd(1);
160    h21->SetFillColor(42);
161    h21->SetXTitle("x");
162    h21->SetYTitle("y");
163    h21->Draw();
164
165    c2->cd(2);
166    hmult->SetFillColor(42);
167    hmult->SetXTitle("multiplicity");
168    hmult->Draw();
169
170    c2->cd(3);
171    hchi1->SetFillColor(42);
172    hchi1->SetXTitle("chi2");
173    hchi1->Draw();
174
175    c2->cd(4);
176    hchi2->SetFillColor(42);
177    hchi2->SetXTitle("chi2");
178    hchi2->Draw();
179
180 }
181
182
183