Merging the VirtualMC branch to the main development branch (HEAD)
[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("galice.root");
13     if (!file) {
14         printf("\n Creating galice.root \n");
15         file = new TFile("galice.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,5.e3);
76    hEzem -> SetXTitle("E (GeV)");
77    TH1F *hPMzem  = new TH1F("hPMzem","Light produced in ZEM PM",100,0,5.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 // NB -> Il TClonesArray delle particelle va inizializzato prima del loop
115    TParticle   *particle;
116    TClonesArray *Particles = gAlice->Particles();
117
118    for (Int_t evNumber=0; evNumber<evTot; evNumber++){
119
120 // Import the Kine and Hits Trees for the event evNumber in the file
121    Int_t nparticles = gAlice->GetEvent(evNumber);
122    if (nparticles <= 0) return;
123    printf("\n   --- nparticles = %d\n",nparticles);
124    
125    Float_t energy, EtotZN=0, EtotZP=0, LightCzn=0, LightCzp=0, LtotZN=0, LtotZP=0;
126    Int_t nbytes=0, nbytesd=0, ipart, nhits, ndigits, pdgcode, ADCzn, ADCzp;
127    AliZDCHit   *ZDChit;
128    AliZDCDigit *ZDCdigit;
129
130 // Get pointers to Alice detectors and Hits containers
131    AliDetector *ZDC  = gAlice->GetDetector("ZDC");
132    if (ZDC) {
133      TClonesArray *ZDChits    = ZDC->Hits();
134 //     TClonesArray *ZDCdigits  = ZDC->Digits();
135    }
136
137 // # of entries in Hits tree
138    TTree *TH = gAlice->TreeH();
139    Int_t ntracks = TH->GetEntries();
140
141 // # of entries in Digits tree
142 //   TTree *TD = gAlice->TreeD();
143 //   Int_t ndigen = TD->GetEntries();
144
145 //   gAlice->ResetDigits();
146 //   nbytesd += TD->GetEvent(ndigen-1);
147
148 //   if (ZDC) {
149 //     ndigits = ZDCdigits->GetEntries();
150 //     printf("\n       Digits Tree --- # of entries: %d; # of digits: %d\n",ndigen, ndigits);
151 //   }
152 //   for(Int_t digit=0; digit<ndigits; digit++) {
153 //     ZDCdigit = (AliZDCDigit*)ZDCdigits->UncheckedAt(digit);
154 //     printf("\n       Digit# %d, fDetector = %d, fVolume = %d, fADCValue = %f\n",
155 //            digit,ZDCdigit->fDetector,ZDCdigit->fQuadrant,ZDCdigit->fADCValue);
156 //   }
157    
158 // Start loop on tracks in the hits containers
159       for (Int_t track=0; track<ntracks; track++) {
160          gAlice->ResetHits();
161          nbytes += TH->GetEvent(track);
162
163          if (ZDC) {
164 //          nhits = ZDChits->GetEntries();
165 //          nhits = ZDChits->GetLast()+1;
166            nhits   = ZDChits->GetEntriesFast();
167 //           printf("\n Hits Tree --- Event %d track %d nhits %d\n",evNumber,track,nhits);
168
169            particle = (TParticle*)Particles->UncheckedAt(track);
170 //         pdgcode = particle->GetPdgCode();
171 //         printf("\nParticle %d\n",pdgcode);
172
173            for(Int_t hit=0; hit<nhits; hit++) {
174              ZDChit   = (AliZDCHit*)ZDChits->UncheckedAt(hit);
175             
176              // Print of the hits
177 //           printf("\nHit # %d, fVolume = %d %d\n",hit,ZDChit->fVolume[0],ZDChit->fVolume[1]);
178 //           printf("Primary energy = %f, Secondary Flag = %d\n",ZDChit->fPrimKinEn,ZDChit->fSFlag);
179 //           printf("Impact point -> %f %f\n",ZDChit->fXImpact,ZDChit->fYImpact);
180 //           printf("Energy = %f, Light in quadrant = %f, Light in common PM = %f\n\n",
181 //                   ZDChit->fEnergy,ZDChit->fLightPMQ,ZDChit->fLightPMC);
182             
183              // Filling histos
184              if(ZDChit->fVolume[0]==1) { //ZN
185                if(ZDChit->fVolume[1]==1)hPMQ1zn->Fill(ZDChit->fLightPMQ);
186                if(ZDChit->fVolume[1]==2)hPMQ2zn->Fill(ZDChit->fLightPMQ);
187                if(ZDChit->fVolume[1]==3)hPMQ3zn->Fill(ZDChit->fLightPMQ);
188                if(ZDChit->fVolume[1]==4)hPMQ4zn->Fill(ZDChit->fLightPMQ);
189                EtotZN   += ZDChit->fEnergy;
190                LtotZN   += (ZDChit->fLightPMQ) + (ZDChit->fLightPMC);
191                LightCzn += ZDChit->fLightPMC;
192                hspotzn->Fill(ZDChit->fXImpact,ZDChit->fYImpact);
193              }
194              if(ZDChit->fVolume[0]==2) { //ZP
195                if(ZDChit->fVolume[1]==1)hPMQ1zp->Fill(ZDChit->fLightPMQ);
196                if(ZDChit->fVolume[1]==2)hPMQ2zp->Fill(ZDChit->fLightPMQ);
197                if(ZDChit->fVolume[1]==3)hPMQ3zp->Fill(ZDChit->fLightPMQ);
198                if(ZDChit->fVolume[1]==4)hPMQ4zp->Fill(ZDChit->fLightPMQ);
199                EtotZP   += ZDChit->fEnergy;
200                LtotZP   += (ZDChit->fLightPMQ) + (ZDChit->fLightPMC);
201                LightCzp += ZDChit->fLightPMC;
202                hspotzp->Fill(ZDChit->fXImpact,ZDChit->fYImpact);
203              }
204              if(ZDChit->fVolume[0]==3) { //ZEM
205                hEzem->Fill(ZDChit->fEnergy);
206                hPMzem->Fill(ZDChit->fLightPMC);
207                hspotzem->Fill(ZDChit->fXImpact,ZDChit->fYImpact);
208              }
209                  
210            }
211            if(nhits!=0){
212              hEzn->Fill(EtotZN);
213              hLzn->Fill(LtotZN);
214              hPMCzn->Fill(LightCzn);
215              hEzp->Fill(EtotZP);
216              hLzp->Fill(LtotZP);
217              hPMCzp->Fill(LightCzp);
218 //             printf("\n Histos var -> Ezn = %f, Lzn = %f, Ezp = %f, Lzp = %f \n\n",
219 //                     EtotZN, LtotZN, EtotZP, LtotZP);
220            }
221          }//ZDC        
222       }//Track loop
223    }//Hit loop 
224    
225 // Control prints
226
227
228 if(detector == 1){ // ZN histos      
229    TCanvas *c1 = new TCanvas("c1","ZN hits",0,10,580,700);
230    c1->cd();
231    TPad *pad1= new TPad("pad1"," ",0.01,0.51,0.99,0.99);
232    TPad *pad2= new TPad("pad2"," ",0.01,0.01,0.49,0.49);
233    TPad *pad3= new TPad("pad3"," ",0.51,0.01,0.99,0.49);
234    pad1->SetFillColor(18);
235    pad2->SetFillColor(18);
236    pad3->SetFillColor(18);
237    pad1->Draw();
238    pad2->Draw();
239    pad3->Draw();
240    pad1->cd();
241    hspotzn->Draw();
242    pad2->cd();
243    hEzn->Draw();
244    pad3->cd();
245    hPMCzn->Draw();
246    
247    TCanvas *c2 = new TCanvas("c2","ZN hits",600,10,600,700);
248    c2->cd();
249    TPad *pad4= new TPad("pad4"," ",0.01,0.51,0.49,0.99);
250    TPad *pad5= new TPad("pad5"," ",0.51,0.51,0.99,0.99);
251    TPad *pad6= new TPad("pad6"," ",0.01,0.01,0.49,0.49);
252    TPad *pad7= new TPad("pad7"," ",0.51,0.01,0.99,0.49);
253    pad4->SetFillColor(18);
254    pad5->SetFillColor(18);
255    pad6->SetFillColor(18);
256    pad7->SetFillColor(18);
257    pad4->Draw();
258    pad5->Draw();
259    pad6->Draw();
260    pad7->Draw();
261    pad4->cd();
262    hPMQ1zn->Draw();
263    pad5->cd();
264    hPMQ2zn->Draw();
265    pad6->cd();
266    hPMQ3zn->Draw();
267    pad7->cd();
268    hPMQ4zn->Draw();
269
270    TCanvas *c3 = new TCanvas("c3","ZN hits",300,10,600,700);
271    c3->cd();
272    hLzn->Draw();
273 }
274
275 if(detector == 2){ // ZP histos      
276    TCanvas *c1 = new TCanvas("c1","ZP hits",0,10,580,700);
277    c1->cd();
278    TPad *pad1= new TPad("pad1"," ",0.01,0.51,0.99,0.99);
279    TPad *pad2= new TPad("pad2"," ",0.01,0.01,0.49,0.49);
280    TPad *pad3= new TPad("pad3"," ",0.51,0.01,0.99,0.49);
281    pad1->SetFillColor(18);
282    pad2->SetFillColor(18);
283    pad3->SetFillColor(18);
284    pad1->Draw();
285    pad2->Draw();
286    pad3->Draw();
287    pad1->cd();
288    hspotzp->Draw();
289    pad2->cd();
290    hEzp->Draw();
291    pad3->cd();
292    hPMCzp->Draw();
293    
294    TCanvas *c2 = new TCanvas("c2","ZP hits",600,10,600,700);
295    c2->cd();
296    TPad *pad4= new TPad("pad4"," ",0.01,0.51,0.49,0.99);
297    TPad *pad5= new TPad("pad5"," ",0.51,0.51,0.99,0.99);
298    TPad *pad6= new TPad("pad6"," ",0.01,0.01,0.49,0.49);
299    TPad *pad7= new TPad("pad7"," ",0.51,0.01,0.99,0.49);
300    pad4->SetFillColor(18);
301    pad5->SetFillColor(18);
302    pad6->SetFillColor(18);
303    pad7->SetFillColor(18);
304    pad4->Draw();
305    pad5->Draw();
306    pad6->Draw();
307    pad7->Draw();
308    pad4->cd();
309    hPMQ1zp->Draw();
310    pad5->cd();
311    hPMQ2zp->Draw();
312    pad6->cd();
313    hPMQ3zp->Draw();
314    pad7->cd();
315    hPMQ4zp->Draw();
316
317    TCanvas *c3 = new TCanvas("c3","ZP hits",300,10,600,700);
318    c3->cd();
319    hLzp->Draw();
320 }
321
322 if(detector == 3){ // ZEM histos      
323    TCanvas *c1 = new TCanvas("c1","ZEM hits",0,10,580,700);
324    c1->cd();
325    TPad *pad1= new TPad("pad1"," ",0.01,0.51,0.99,0.99);
326    TPad *pad2= new TPad("pad2"," ",0.01,0.01,0.99,0.49);
327    pad1->SetFillColor(18);
328    pad2->SetFillColor(18);
329    pad1->Draw();
330    pad2->Draw();
331    pad1->cd();
332    hspotzem->Draw();
333    pad2->cd();
334    hEzem->Draw();
335    
336    TCanvas *c2 = new TCanvas("c2","ZEM hits",600,10,600,700);
337    c2->cd();
338    hPMzem->Draw();
339 }
340 //   file->Close();   
341 }