Macro to read Hits and Digits Trees
[u/mrichter/AliRoot.git] / ZDC / ZDCtest.C
1 void ZDCtest (Int_t detector=0, Int_t evTot = 0) 
2 {
3    delete gAlice;
4    gAlice=0;
5 // Dynamically link some shared libs
6    if (gClassTable->GetID("AliRun") < 0) {
7       gROOT->LoadMacro("loadlibs.C");
8       loadlibs();
9    }
10       
11 // Connect the Root Galice file containing Geometry, Kine, Hits and Digits
12    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("Hijing_b2.root");
13     if (!file) {
14         printf("\n Creating galice.root \n");
15         file = new TFile("Hijing_b2.root");
16     } else {
17         printf("\n galice.root found in file list");
18     }
19     file->ls();
20
21 // Get AliRun object from file or create it if not on file
22    if (!gAlice) {
23       gAlice = (AliRun*)file->Get("gAlice");
24       if (gAlice) printf("AliRun object found on file\n");
25       if (!gAlice) {
26             printf("\n create new gAlice object");
27             gAlice = new AliRun("gAlice","Alice test program");
28         }
29    }
30    file->ls();
31    
32 // Create some histograms
33
34    TH2F *hspotzn  = new TH2F("hspotzn","Y vs X on ZN front face",100,-4.,4.,100,-4.,4.);
35    hspotzn -> SetXTitle("X on ZN");
36    hspotzn -> SetYTitle("Y on ZN");
37    TH2F *hspotzp  = new TH2F("hspotzp","Y vs X on ZP front face",100,-12.,12.,100,-12.,12.);
38    hspotzp -> SetXTitle("X on ZP");
39    hspotzp -> SetYTitle("Y on ZP");
40    TH2F *hspotzem = new TH2F("hspotzem","Y vs X on ZEM front face",100,-4.,4.,100,-4.,4.);
41    hspotzem -> SetXTitle("X on ZEM");
42    hspotzem -> SetYTitle("Y on ZEM");
43
44    TH1F *hEzn    = new TH1F("hEzn","Energy deposited in ZN",100,0,4.e3);
45    hEzn -> SetXTitle("E (GeV)");
46    TH1F *hLzn    = new TH1F("hLzn", "Total light in ZN    ",100,0,4.e3);
47    hLzn -> SetXTitle("phe");
48    TH1F *hPMCzn  = new TH1F("hPMCzn", "Light in common PM    ",100,0,1.e3);
49    hPMCzn -> SetXTitle("phe");
50    TH1F *hPMQ1zn = new TH1F("hPMQ1zn","Light in quadrant 1 PM",100,0,1.e3);
51    hPMQ1zn -> SetXTitle("phe");
52    TH1F *hPMQ2zn = new TH1F("hPMQ2zn","Light in quadrant 2 PM",100,0,1.e3);
53    hPMQ2zn -> SetXTitle("phe");
54    TH1F *hPMQ3zn = new TH1F("hPMQ3zn","Light in quadrant 3 PM",100,0,1.e3);
55    hPMQ3zn -> SetXTitle("phe");
56    TH1F *hPMQ4zn = new TH1F("hPMQ4zn","Light in quadrant 4 PM",100,0,1.e3);
57    hPMQ4zn -> SetXTitle("phe");
58
59    TH1F *hEzp    = new TH1F("hEzp","Energy deposited in ZP",100,0,4.e3);
60    hEzp -> SetXTitle("E (GeV)");
61    TH1F *hLzp    = new TH1F("hLzp", "Total light in ZP    ",100,0,4.e3);
62    hLzp -> SetXTitle("phe");
63    TH1F *hPMCzp  = new TH1F("hPMCzp","Light in common PM     ",100,0,1.e3);
64    hPMCzp -> SetXTitle("phe");
65    TH1F *hPMQ1zp = new TH1F("hPMQ1zp","Light in quadrant 1 PM",100,0,1.e3);
66    hPMQ1zp -> SetXTitle("phe");
67    TH1F *hPMQ2zp = new TH1F("hPMQ2zp","Light in quadrant 2 PM",100,0,1.e3);
68    hPMQ2zp -> SetXTitle("phe");
69    TH1F *hPMQ3zp = new TH1F("hPMQ3zp","Light in quadrant 3 PM",100,0,1.e3);
70    hPMQ3zp -> SetXTitle("phe");
71    TH1F *hPMQ4zp = new TH1F("hPMQ4zp","Light in quadrant 4 PM",100,0,1.e3);
72    hPMQ4zp -> SetXTitle("phe");
73
74
75    TH1F *hEzem   = new TH1F("hEzem","Energy deposited in ZEM",100,0,1.e3);
76    hEzem -> SetXTitle("E (GeV)");
77    TH1F *hPMzem  = new TH1F("hPMzem","Light produced in ZEM PM",100,0,4.e3);
78    hPMzem -> SetXTitle("phe");
79    
80    //
81    
82    TH1F *dPMCzn  = new TH1F("dPMCzn","Common PM     ",100,0,1.e3);
83    dPMCzn -> SetXTitle("ADC channels");
84    TH1F *dPMQ1zn = new TH1F("dPMQ1zn","Quadrant 1 PM",100,0,1.e3);
85    dPMQ1zn -> SetXTitle("ADC channels");
86    TH1F *dPMQ2zn = new TH1F("dPMQ2zn","Quadrant 2 PM",100,0,1.e3);
87    dPMQ2zn -> SetXTitle("ADC channels");
88    TH1F *dPMQ3zn = new TH1F("dPMQ3zn","Quadrant 3 PM",100,0,1.e3);
89    dPMQ3zn -> SetXTitle("ADC channels");
90    TH1F *dPMQ4zn = new TH1F("dPMQ4zn","Quadrant 4 PM",100,0,1.e3);
91    dPMQ4zn -> SetXTitle("ADC channels");
92    TH1F *dZN     = new TH1F("dZN","Total light in ZN",100,0,1.e3);
93    dZN -> SetXTitle("ADC channels");
94    
95    TH1F *dPMCzp  = new TH1F("dPMCzp","Common PM     ",100,0,1.e3);
96    dPMCzp -> SetXTitle("ADC channels");
97    TH1F *dPMQ1zp = new TH1F("dPMQ1zp","Quadrant 1 PM",100,0,1.e3);
98    dPMQ1zp -> SetXTitle("ADC channels");
99    TH1F *dPMQ2zp = new TH1F("dPMQ2zp","Quadrant 2 PM",100,0,1.e3);
100    dPMQ2zp -> SetXTitle("ADC channels");
101    TH1F *dPMQ3zp = new TH1F("dPMQ3zp","Quadrant 3 PM",100,0,1.e3);
102    dPMQ3zp -> SetXTitle("ADC channels");
103    TH1F *dPMQ4zp = new TH1F("dPMQ4zp","Quadrant 4 PM",100,0,1.e3);
104    dPMQ4zp -> SetXTitle("ADC channels");
105    TH1F *dZP     = new TH1F("dZP","Total light in ZP",100,0,1.e3);
106    dZP -> SetXTitle("ADC channels");
107
108    TH1F *dZEM     = new TH1F("dZEM","Total light in ZEM",100,0,1.e3);
109    dZEM -> SetXTitle("ADC channels");
110
111 //
112 //   Loop over events 
113 //
114
115    for (Int_t evNumber=0; evNumber<evTot; evNumber++){
116
117 // Import the Kine and Hits Trees for the event evNumber in the file
118    Int_t nparticles = gAlice->GetEvent(evNumber);
119    if (nparticles <= 0) return;
120    printf("\n   --- nparticles = %d\n",nparticles);
121    
122    Float_t energy, EtotZN=0, EtotZP=0, LightCzn=0, LightCzp=0, LtotZN=0, LtotZP=0;
123    Int_t nbytes=0, nbytesd=0, ipart, nhits, ndigits, pdgcode, ADCzn, ADCzp;
124    TParticle   *particle;
125    AliZDCHit   *ZDChit;
126    AliZDCDigit *ZDCdigit;
127
128 // Get pointers to Alice detectors and Hits containers
129    AliDetector *ZDC  = gAlice->GetDetector("ZDC");
130    TClonesArray *Particles = gAlice->Particles();
131    if (ZDC) {
132      TClonesArray *ZDChits    = ZDC->Hits();
133      TClonesArray *ZDCdigits  = ZDC->Digits();
134    }
135
136 // # of entries in Hits tree
137    TTree *TH = gAlice->TreeH();
138    Int_t ntracks = TH->GetEntries();
139
140 // # of entries in Digits tree
141    TTree *TD = gAlice->TreeD();
142    Int_t ndigen = TD->GetEntries();
143
144    gAlice->ResetDigits();
145    nbytesd += TD->GetEvent(ndigen-1);
146
147    if (ZDC) {
148      ndigits = ZDCdigits->GetEntries();
149      printf("\n Digits Tree --- # of entries: %d; # of digits: %d\n",ndigen, ndigits);
150    }
151    for(Int_t digit=0; digit<ndigits; digit++) {
152      ZDCdigit = (AliZDCDigit*)ZDCdigits->UncheckedAt(digit);
153      printf("\n Digit# %d, fDetector = %d, fVolume = %d, fADCValue = %f\n",
154             digit,ZDCdigit->fDetector,ZDCdigit->fQuadrant,ZDCdigit->fADCValue);
155    }
156    
157 // Start loop on tracks in the hits containers
158       for (Int_t track=0; track<ntracks; track++) {
159          gAlice->ResetHits();
160          nbytes += TH->GetEvent(track);
161
162          if (ZDC) {
163 //          nhits = ZDChits->GetEntries();
164 //          nhits = ZDChits->GetLast()+1;
165            nhits   = ZDChits->GetEntriesFast();
166 //           printf("\n Hits Tree --- Event %d track %d nhits %d\n",evNumber,track,nhits);
167
168            particle = (TParticle*)Particles->UncheckedAt(track);
169 //         pdgcode = particle->GetPdgCode();
170 //         printf("\nParticle %d\n",pdgcode);
171
172            for(Int_t hit=0; hit<nhits; hit++) {
173              ZDChit   = (AliZDCHit*)ZDChits->UncheckedAt(hit);
174             
175              // Print of the hits
176 //           printf("\nHit # %d, fVolume = %d %d\n",hit,ZDChit->fVolume[0],ZDChit->fVolume[1]);
177 //           printf("Primary energy = %f, Secondary Flag = %d\n",ZDChit->fPrimKinEn,ZDChit->fSFlag);
178 //           printf("Impact point -> %f %f\n",ZDChit->fXImpact,ZDChit->fYImpact);
179 //           printf("Energy = %f, Light in quadrant = %f, Light in common PM = %f\n\n",
180 //                   ZDChit->fEnergy,ZDChit->fLightPMQ,ZDChit->fLightPMC);
181             
182              // Filling histos
183              if(ZDChit->fVolume[0]==1) { //ZN
184                if(ZDChit->fVolume[1]==1)hPMQ1zn->Fill(ZDChit->fLightPMQ);
185                if(ZDChit->fVolume[1]==2)hPMQ2zn->Fill(ZDChit->fLightPMQ);
186                if(ZDChit->fVolume[1]==3)hPMQ3zn->Fill(ZDChit->fLightPMQ);
187                if(ZDChit->fVolume[1]==4)hPMQ4zn->Fill(ZDChit->fLightPMQ);
188                EtotZN   += ZDChit->fEnergy;
189                LtotZN   += (ZDChit->fLightPMQ) + (ZDChit->fLightPMC);
190                LightCzn += ZDChit->fLightPMC;
191                hspotzn->Fill(ZDChit->fXImpact,ZDChit->fYImpact);
192              }
193              if(ZDChit->fVolume[0]==2) { //ZP
194                if(ZDChit->fVolume[1]==1)hPMQ1zp->Fill(ZDChit->fLightPMQ);
195                if(ZDChit->fVolume[1]==2)hPMQ2zp->Fill(ZDChit->fLightPMQ);
196                if(ZDChit->fVolume[1]==3)hPMQ3zp->Fill(ZDChit->fLightPMQ);
197                if(ZDChit->fVolume[1]==4)hPMQ4zp->Fill(ZDChit->fLightPMQ);
198                EtotZP   += ZDChit->fEnergy;
199                LtotZP   += (ZDChit->fLightPMQ) + (ZDChit->fLightPMC);
200                LightCzp += ZDChit->fLightPMC;
201                hspotzp->Fill(ZDChit->fXImpact,ZDChit->fYImpact);
202              }
203              if(ZDChit->fVolume[0]==3) { //ZEM
204                hEzem->Fill(ZDChit->fEnergy);
205                hPMzem->Fill(ZDChit->fLightPMQ);
206                hspotzem->Fill(ZDChit->fXImpact,ZDChit->fYImpact);
207              }
208                  
209            }
210            if(nhits!=0){
211              hEzn->Fill(EtotZN);
212              hLzn->Fill(LtotZN);
213              hPMCzn->Fill(LightCzn);
214              hEzp->Fill(EtotZP);
215              hLzp->Fill(LtotZP);
216              hPMCzp->Fill(LightCzp);
217              printf("\n Histos var -> Ezn = %f, Lzn = %f, Ezp = %f, Lzp = %f \n\n",
218                      EtotZN, LtotZN, EtotZP, LtotZP);
219            }
220          }//ZDC        
221       }//Track loop
222    }//Hit loop 
223    
224 // Control prints
225
226
227 if(detector == 1){ // ZN histos      
228    TCanvas *c1 = new TCanvas("c1","ZN hits",0,10,580,700);
229    c1->cd();
230    TPad *pad1= new TPad("pad1"," ",0.01,0.51,0.99,0.99);
231    TPad *pad2= new TPad("pad2"," ",0.01,0.01,0.49,0.49);
232    TPad *pad3= new TPad("pad3"," ",0.51,0.01,0.99,0.49);
233    pad1->SetFillColor(18);
234    pad2->SetFillColor(18);
235    pad3->SetFillColor(18);
236    pad1->Draw();
237    pad2->Draw();
238    pad3->Draw();
239    pad1->cd();
240    hspotzn->Draw();
241    pad2->cd();
242    hEzn->Draw();
243    pad3->cd();
244    hPMCzn->Draw();
245    
246    TCanvas *c2 = new TCanvas("c2","ZN hits",600,10,600,700);
247    c2->cd();
248    TPad *pad4= new TPad("pad4"," ",0.01,0.51,0.49,0.99);
249    TPad *pad5= new TPad("pad5"," ",0.51,0.51,0.99,0.99);
250    TPad *pad6= new TPad("pad6"," ",0.01,0.01,0.49,0.49);
251    TPad *pad7= new TPad("pad7"," ",0.51,0.01,0.99,0.49);
252    pad4->SetFillColor(18);
253    pad5->SetFillColor(18);
254    pad6->SetFillColor(18);
255    pad7->SetFillColor(18);
256    pad4->Draw();
257    pad5->Draw();
258    pad6->Draw();
259    pad7->Draw();
260    pad4->cd();
261    hPMQ1zn->Draw();
262    pad5->cd();
263    hPMQ2zn->Draw();
264    pad6->cd();
265    hPMQ3zn->Draw();
266    pad7->cd();
267    hPMQ4zn->Draw();
268
269    TCanvas *c3 = new TCanvas("c3","ZN hits",300,10,600,700);
270    c3->cd();
271    hLzn->Draw();
272 }
273
274 if(detector == 2){ // ZP histos      
275    TCanvas *c1 = new TCanvas("c1","ZP hits",0,10,580,700);
276    c1->cd();
277    TPad *pad1= new TPad("pad1"," ",0.01,0.51,0.99,0.99);
278    TPad *pad2= new TPad("pad2"," ",0.01,0.01,0.49,0.49);
279    TPad *pad3= new TPad("pad3"," ",0.51,0.01,0.99,0.49);
280    pad1->SetFillColor(18);
281    pad2->SetFillColor(18);
282    pad3->SetFillColor(18);
283    pad1->Draw();
284    pad2->Draw();
285    pad3->Draw();
286    pad1->cd();
287    hspotzp->Draw();
288    pad2->cd();
289    hEzp->Draw();
290    pad3->cd();
291    hPMCzp->Draw();
292    
293    TCanvas *c2 = new TCanvas("c2","ZP hits",600,10,600,700);
294    c2->cd();
295    TPad *pad4= new TPad("pad4"," ",0.01,0.51,0.49,0.99);
296    TPad *pad5= new TPad("pad5"," ",0.51,0.51,0.99,0.99);
297    TPad *pad6= new TPad("pad6"," ",0.01,0.01,0.49,0.49);
298    TPad *pad7= new TPad("pad7"," ",0.51,0.01,0.99,0.49);
299    pad4->SetFillColor(18);
300    pad5->SetFillColor(18);
301    pad6->SetFillColor(18);
302    pad7->SetFillColor(18);
303    pad4->Draw();
304    pad5->Draw();
305    pad6->Draw();
306    pad7->Draw();
307    pad4->cd();
308    hPMQ1zp->Draw();
309    pad5->cd();
310    hPMQ2zp->Draw();
311    pad6->cd();
312    hPMQ3zp->Draw();
313    pad7->cd();
314    hPMQ4zp->Draw();
315
316    TCanvas *c3 = new TCanvas("c3","ZP hits",300,10,600,700);
317    c3->cd();
318    hLzp->Draw();
319 }
320
321 if(detector == 3){ // ZEM histos      
322    TCanvas *c1 = new TCanvas("c1","ZEM hits",0,10,580,700);
323    c1->cd();
324    TPad *pad1= new TPad("pad1"," ",0.01,0.51,0.99,0.99);
325    TPad *pad2= new TPad("pad2"," ",0.01,0.01,0.99,0.49);
326    pad1->SetFillColor(18);
327    pad2->SetFillColor(18);
328    pad1->Draw();
329    pad2->Draw();
330    pad1->cd();
331    hspotzem->Draw();
332    pad2->cd();
333    hEzem->Draw();
334    
335    TCanvas *c2 = new TCanvas("c2","ZEM hits",600,10,600,700);
336    c2->cd();
337    hPMzem->Draw();
338 }
339 //   file->Close();   
340 }