New TTask based method to do Digits To clusters. Works with files of multiple
[u/mrichter/AliRoot.git] / MUON / MUONocc.C
1 void MUONocc (Int_t evNumber1=0,Int_t evNumber2=0, Int_t ic=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 // Dynamically link some shared libs                    
11
12    if (gClassTable->GetID("AliRun") < 0) {
13       gROOT->LoadMacro("loadlibs.C");
14       loadlibs();
15    }
16 // Connect the Root Galice file containing Geometry, Kine and Hits
17
18     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
19
20     
21     if (!file) {
22         printf("\n Creating galice.root \n");
23         file = new TFile("galice.root");
24     } else {
25         printf("\n galice.root found in file list");
26     }
27     file->ls();
28     
29 // Get AliRun object from file or create it if not on file
30     if (!gAlice) {
31         gAlice = (AliRun*)(file->Get("gAlice"));
32         if (gAlice) printf("AliRun object found on file\n");
33         if (!gAlice) {
34             printf("\n create new gAlice object");
35             gAlice = new AliRun("gAlice","Alice test program");
36         }
37     }
38     
39 //  Create some histograms
40
41     TH1F *Hits[10], *HitDensity[10], *Occ0[10],*Occ1[10], *Mult[10];
42     for (Int_t i=0; i<10; i++) {
43         Hits[i]     =  new TH1F("hits1","Hit Distribution",26,0,260);
44         Hits[i]->SetFillColor(0);
45         Hits[i]->SetXTitle("R (cm)");
46
47         HitDensity[i]  =  new TH1F("dhits1","Hit Density Distribution",26,0,260);
48         HitDensity[i]->SetFillColor(0);
49         HitDensity[i]->SetXTitle("R (cm)");
50
51         Occ0[i]  =  new TH1F("occ0","Occupancy Density",26,0,260);
52         Occ0[i] -> SetFillColor(0);
53         Occ0[i] -> SetXTitle("R (cm)");
54         Occ1[i]  =  new TH1F("occ1","Occupancy",26,0,260);
55         Occ1[i] -> SetFillColor(0);
56         Occ1[i] -> SetXTitle("R (cm)");
57         Occ1[i] -> SetYTitle("Occupancy");
58         Mult[i]  =  new TH1F("mult","Mult distribution",26,0,260);
59         Mult[i] -> SetFillColor(0);
60         Mult[i] -> SetXTitle("R (cm)");
61     }
62
63
64    TH1F *theta   =  new TH1F("theta","Theta distribution",180,0,180);
65  
66    TH1F *emult   =  new TH1F("emult","Event Multiplicity",100,0,1000);   
67    AliMUONChamber*  iChamber;
68    AliMUONSegmentation* seg;
69    
70 //
71 //   Loop over events 
72 //
73    Float_t dpx[2], dpy[2];
74    Int_t mdig =0;
75    Int_t nhit=0;
76    
77    for (Int_t nev=0; nev<= evNumber2; nev++) {
78        Int_t nparticles = gAlice->GetEvent(nev);
79        if (nev < evNumber1) continue;
80        if (nparticles <= 0) return;
81
82        AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
83        printf("\n track %d %d \n ", nev, MUON);
84
85        TTree *TH = gAlice->TreeH();
86        Int_t ntracks = TH->GetEntries();
87 //
88 //   Loop over events
89 //
90        Int_t EvMult=0;
91        
92        for (Int_t track=0; track<ntracks;track++) {
93            gAlice->ResetHits();
94            Int_t nbytes += TH->GetEvent(track);
95            if (MUON)  {
96 //
97 //  Loop over hits
98 //
99                Int_t NperCh=0;
100                
101                for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1); 
102                    mHit;
103                    mHit=(AliMUONHit*)MUON->NextHit()) 
104                {
105                    Int_t   nch   = mHit->Chamber();  // chamber number
106                    Float_t x     = mHit->X();        // x-pos of hit
107                    Float_t y     = mHit->Y();        // y-pos
108                    Float_t Eloss = mHit->Eloss();
109                    Float_t Theta = mHit->Theta();
110                    theta->Fill(Theta,(float) 1);
111                    if (nch >10) continue;
112                    if (nch ==1) EvMult++;
113 //
114 //
115                    iChamber = & MUON->Chamber(nch-1);
116                    response= iChamber.GetResponseModel();
117                    seg     = iChamber.GetSegmentationModel(1);
118                    NperCh++;
119                    Int_t i,j;
120                    seg->GetPadIxy(x,y,i,j);
121                    Int_t isec = seg->Sector(i,j);
122                    if (isec ==1 && nch==ic) nhit++;
123                    
124                    Float_t a=seg->Dpx(isec)*seg->Dpy(isec);
125                    Float_t r=TMath::Sqrt(x*x+y*y);
126                    Float_t wgt=1/(2*10*TMath::Pi()*r)/(evNumber2+1);
127                    wgt=wgt*(1.+24./(2.*TMath::Pi()*r));
128                    Hits[nch-1]->Fill(r,(float) 1);
129                    HitDensity[nch-1]->Fill(r,wgt);
130                    Occ0[nch-1]->Fill(r,wgt*a);
131                } // hit loop
132            } // if MUON
133        } // track loop
134        emult->Fill(Float_t(EvMult), (Float_t) 1);
135        
136        Int_t iseg=1;
137
138        for (Int_t ich=0; ich<10; ich++) {
139            iChamber = & MUON->Chamber(ich);
140            seg=iChamber.GetSegmentationModel(iseg);
141            gAlice->ResetDigits();
142            gAlice->TreeD()->GetEvent(iseg);
143            TClonesArray *MUONDigits  = MUON->DigitsAddress(ich);
144            Int_t Ndigits=MUONDigits->GetEntriesFast();
145            AliMUONDigit* dig;
146            printf("\n Reading %d digits\n", Ndigits);
147            if (MUONDigits)  {
148                for (Int_t ndig=0; ndig<Ndigits; ndig++) 
149                {
150                    dig = (AliMUONDigit*)MUONDigits->UncheckedAt(ndig); 
151                    Int_t i=dig->PadX();
152                    Int_t j=dig->PadY();
153                    Float_t x,y;
154                    seg->GetPadCxy(i,j,x,y);
155                    Int_t isec = seg->Sector(i,j);
156                    Float_t a=seg->Dpx(isec)*seg->Dpy(isec);
157                    Float_t r=TMath::Sqrt(x*x+y*y);
158 //             if (r<25) 
159 //             printf("\n Sector,a %d %f", isec,a);
160                    if (isec==1) mdig++;
161                    
162                    Float_t wgt;
163                    wgt=1/(2.*10.*TMath::Pi()*r)/(evNumber2+1)*a;
164 // Take into account inefficiency due to frames
165                    wgt=wgt*(1.+24./(2.*TMath::Pi()*r));
166                    Occ1[ich]->Fill(r,wgt);
167                } // digit loop
168                Mult[ich]->Divide(Occ1[ich],Occ0[ich]);
169            } // chamber loop
170        } // if MUONDigits
171    } // event loop 
172 //
173    printf("\n hits, digits %d %d\n ", nhit, mdig);
174    
175
176    
177 //Create a canvas, set the view range, show histograms
178    TCanvas *c1 = new TCanvas("c1","Hit Densities",400,10,600,700);
179    pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
180    pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
181    pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
182    pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
183    pad11->SetFillColor(0);
184    pad12->SetFillColor(0);
185    pad13->SetFillColor(0);
186    pad14->SetFillColor(0);
187    pad11->Draw();
188    pad12->Draw();
189    pad13->Draw();
190    pad14->Draw();
191    
192    pad11->cd();
193    HitDensity[0]->Draw();
194    pad12->cd();
195    HitDensity[1]->Draw();
196    pad13->cd();
197    HitDensity[2]->Draw();
198    pad14->cd();
199    HitDensity[3]->Draw();
200
201    TCanvas *c2 = new TCanvas("c2","Hit Densities",400,10,600,700);
202    pad21 = new TPad("pad21"," ",0.01,0.51,0.49,0.99);
203    pad22 = new TPad("pad22"," ",0.51,0.51,0.99,0.99);
204    pad23 = new TPad("pad23"," ",0.01,0.01,0.49,0.49);
205    pad24 = new TPad("pad24"," ",0.51,0.01,0.99,0.49);
206    pad21->SetFillColor(0);
207    pad22->SetFillColor(0);
208    pad23->SetFillColor(0);
209    pad24->SetFillColor(0);
210    pad21->Draw();
211    pad22->Draw();
212    pad23->Draw();
213    pad24->Draw();
214
215    pad21->cd();
216    HitDensity[4]->Draw();
217    pad22->cd();
218    HitDensity[5]->Draw();
219    pad23->cd();
220    HitDensity[6]->Draw();
221    pad24->cd();
222    HitDensity[7]->Draw();
223
224    TCanvas *c3 = new TCanvas("c3","Hit Densities",400,10,600,700);
225    pad31 = new TPad("pad31"," ",0.01,0.51,0.49,0.99);
226    pad32 = new TPad("pad32"," ",0.51,0.51,0.99,0.99);
227    pad33 = new TPad("pad33"," ",0.01,0.01,0.49,0.49);
228    pad34 = new TPad("pad34"," ",0.51,0.01,0.99,0.49);
229    pad31->SetFillColor(0);
230    pad32->SetFillColor(0);
231    pad33->SetFillColor(0);
232    pad34->SetFillColor(0);
233    pad31->Draw();
234    pad32->Draw();
235    pad33->Draw();
236    pad34->Draw();
237
238    pad31->cd();
239    HitDensity[8]->Draw();
240    pad32->cd();
241    HitDensity[9]->Draw();
242
243    TCanvas *c4 = new TCanvas("c4","Occupancies",400,10,600,700);
244    pad41 = new TPad("pad41"," ",0.01,0.51,0.49,0.99);
245    pad42 = new TPad("pad42"," ",0.51,0.51,0.99,0.99);
246    pad43 = new TPad("pad43"," ",0.01,0.01,0.49,0.49);
247    pad44 = new TPad("pad44"," ",0.51,0.01,0.99,0.49);
248    pad41->SetFillColor(0);
249    pad42->SetFillColor(0);
250    pad43->SetFillColor(0);
251    pad44->SetFillColor(0);
252    pad41->Draw();
253    pad42->Draw();
254    pad43->Draw();
255    pad44->Draw();
256
257
258    pad41->cd();
259    Occ1[0]->Draw();
260    pad42->cd();
261    Occ1[2]->Draw();
262    pad43->cd();
263    Occ1[4]->Draw();
264    pad44->cd();
265    Occ1[6]->Scale(1.25);
266    Occ1[6]->Draw();
267
268    TCanvas *c5 = new TCanvas("c5","Occupancies",400,10,600,700);
269    pad51 = new TPad("pad41"," ",0.01,0.51,0.49,0.99);
270    pad52 = new TPad("pad42"," ",0.51,0.51,0.99,0.99);
271    pad53 = new TPad("pad43"," ",0.01,0.01,0.49,0.49);
272    pad54 = new TPad("pad44"," ",0.51,0.01,0.99,0.49);
273    pad51->SetFillColor(0);
274    pad52->SetFillColor(0);
275    pad53->SetFillColor(0);
276    pad54->SetFillColor(0);
277    pad51->Draw();
278    pad52->Draw();
279    pad53->Draw();
280    pad54->Draw();
281
282
283    pad51->cd();
284    Occ1[8]->Scale(1.25);
285    
286    Occ1[8]->Draw();
287 }
288
289
290
291
292
293
294
295
296
297