30bbbe40f92d61399ae6a13537dae977a598ab4e
[u/mrichter/AliRoot.git] / macros / analHits.C
1 void analHits (const char *filename="galice.root",Int_t evNumber=0, char *opt="Liny"){
2 /////////////////////////////////////////////////////////////////////////
3 //   This macro is a small example of a ROOT macro
4 //   illustrating how to read the output of GALICE
5 //   and fill some histograms.
6 //   
7 //     Root > .L analHits.C               //this loads the macro in memory
8 //     Root > analHits();                 //by default process first event   
9 //     Root > analHits("galice2.root",2); //process third event from 
10 //                                          galice2.root file.
11 //Begin_Html
12 /*
13 <img src="picts/analHits.gif">
14 */
15 //End_Html
16 /////////////////////////////////////////////////////////////////////////
17   if(gAlice){
18     delete gAlice;
19     gAlice=0;
20   }
21   else{
22     // Dynamically link some shared libs
23     if(gClassTable->GetID("AliRun") < 0) {
24       gROOT->LoadMacro("loadlibs.C");
25       loadlibs();
26     } // end if
27   }
28 // Connect the Root Galice file containing Geometry, Kine and Hits
29     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject(filename);
30     if(!file) file = new TFile(filename);
31
32 // Get AliRun object from file or create it if not on file
33     if(!gAlice) {
34         gAlice = (AliRun*)file->Get("gAlice");
35         if(gAlice) printf("AliRun object found on file\n");
36         if(!gAlice) gAlice = new AliRun("gAlice","Alice test program");
37     } // end if !gAlice
38       
39 // Set event pointer to this event
40     Int_t nparticles = gAlice->GetEvent(evNumber);
41     if (nparticles <= 0){
42         cout << "No particles found for event " << evNumber;
43         cout << " in file " << filename << endl;
44         return;
45     } // end if nparticles <=0
46
47 // Pointer to specific detector hits.
48     AliFMDhit    *fmdHit;
49     AliITShit    *itsHit;
50     AliMUONHit   *muonHit;
51     AliPHOSHit   *phosHit;
52     AliPMDhit    *pmdHit;
53     AliRICHHit   *richHit;
54     AliSTARThit  *startHit;
55     AliTOFhit    *tofHit;
56     AliTPChit    *tpcHit;
57     AliTPCTrackHits *tpc2Hit;
58     AliTRDhit    *trdHit;
59     AliCASTORhit *castorHit;
60     AliZDCHit    *zdcHit;
61
62 // Get pointers to ALL Alice detectors and Hits containers
63     AliFMD    *FMD    = (AliFMD*)    gAlice->GetDetector("FMD");
64     AliITS    *ITS    = (AliITS*)    gAlice->GetDetector("ITS");
65     AliMUON   *MUON   = (AliMUON*)   gAlice->GetDetector("MUON");
66     AliPHOS   *PHOS   = (AliPHOS*)   gAlice->GetDetector("PHOS");
67     AliPMD    *PMD    = (AliPMD*)    gAlice->GetDetector("PMD");
68     AliRICH   *RICH   = (AliRICH*)   gAlice->GetDetector("RICH");
69     AliSTART  *START  = (AliSTART*)  gAlice->GetDetector("START");
70     AliTOF    *TOF    = (AliTOF*)    gAlice->GetDetector("TOF");
71     AliTPC    *TPC    = (AliTPC*)    gAlice->GetDetector("TPC");
72     AliTPC    *TPC    = (AliTPC*)    gAlice->GetDetector("TPC");
73     AliTRD    *TRD    = (AliTRD*)    gAlice->GetDetector("TRD");
74     AliCASTOR *CASTOR = (AliCASTOR*) gAlice->GetDetector("CASTOR");
75     AliZDC    *ZDC    = (AliZDC*)    gAlice->GetDetector("ZDC");
76
77 // Get pointer to the particles
78 //    TClonesArray *Particles = gAlice->Particles();
79 //    TParticle    *part;
80
81     // Create histograms
82     if(FMD)   TH1F *hFMD   = new TH1F("hFMD"   ,"Hit Radius",100,0.,100.);
83     if(ITS)   TH1F *hITS   = new TH1F("hITS"   ,"Ionization",100,0.,3.e-3);
84     if(MUON)  TH1F *hMUON  = new TH1F("hMUON"  ,"Hit Radius",100,0.,500.);
85     if(PHOS)   TH1F *hPHOS  = new TH1F("hPHOS"  ,"Energy Dep.",100,0.,0.5);
86     if(PMD)     TH1F *hPMD   = new TH1F("hPMD"   ,"Energy Dep.",100,0.,1.e+5);
87     if(RICH)  TH1F *hRICH  = new TH1F("hRICH"  ,"Energy loss",100,0.,1.e-5);
88     if(START) TH1F *hSTART = new TH1F("hSTART" ,"Time of Flight",100,0.,10.);
89     if(TOF)   TH1F *hTOF   = new TH1F("hTOF"   ,"Time of Flight",100,0.,1.e-5);
90     if(TPC)   TH1F *hTPC   = new TH1F("hTPC"   ,"Charge",100,0.,70.2);
91     if(TRD)   TH1F *hTRD   = new TH1F("hTRD"   ,"Charge",100,0.,10.);
92     if(CASTOR)TH1F *hCASTOR= new TH1F("hCASTOR","Hit Radius",100,0.,10.);
93     if(ZDC)   TH1F *hZDC   = new TH1F("hZDC"   ,"Energy",100,0.,5.);
94 //    TH1F *hTPAR  = new TH1F("hTPAR" ,"?",6,1,7);  
95     Int_t track,ntracks = gAlice->TreeH()->GetEntries();
96 // Start loop on tracks in the hits containers
97     for(track=0; track<ntracks;track++){
98       //MI change
99       gAlice->ResetHits();
100       gAlice->TreeH()->GetEvent(track);
101       if(FMD){
102             for(fmdHit=(AliFMDhit*)FMD->FirstHit(-1);fmdHit;
103                 fmdHit=(AliFMDhit*)FMD->NextHit()){
104                 hFMD->Fill(TMath::Hypot(fmdHit->X(),fmdHit->Y()));
105             } // end for fmdHit
106         } // end if FMD
107         if(ITS){
108             for(itsHit=(AliITShit*)ITS->FirstHit(-1);itsHit;
109                 itsHit=(AliITShit*)ITS->NextHit()){
110                 if(itsHit->GetIonization()>0.0){//only after a step in the ITS
111                     hITS->Fill(itsHit->GetIonization());
112                 } // end if
113             } // end for itsHit
114         } // end if ITS
115         if(MUON){
116             for(muonHit=(AliMUONHit*)MUON->FirstHit(-1);muonHit;
117                 muonHit=(AliMUONHit*)MUON->NextHit()){
118                 hMUON->Fill(TMath::Hypot(muonHit->X(),muonHit->Y()));
119             } // end for muonHit
120         } // end if MUON
121         if(PHOS){
122             for(phosHit=(AliPHOSHit*)PHOS->FirstHit(-1);phosHit;
123                 phosHit=(AliPHOSHit*)PHOS->NextHit()){
124                 hPHOS->Fill(phosHit->GetEnergy());
125             } // end for phosHit
126         } // end if PHOS
127         if(PMD){
128             for(pmdHit=(AliPMDhit*)PMD->FirstHit(-1);pmdHit;
129                 pmdHit=(AliPMDhit*)PMD->NextHit()){
130                 hPMD->Fill(pmdHit->GetEnergy());
131             } // end for pmdHit
132         } // end if PMD
133         if(RICH){
134             for(richHit=(AliRICHHit*)RICH->FirstHit(-1);richHit;
135                 richHit=(AliRICHHit*)RICH->NextHit()){
136                 hRICH->Fill(richHit->fEloss);
137             } // end for richHit
138         } // end if RICH
139         if(START){
140             for(startHit=(AliSTARThit*)START->FirstHit(-1);startHit;
141                 startHit=(AliSTARThit*)START->NextHit()){
142                 hSTART->Fill(startHit->fTime);
143             } // end for startHit
144         } // end if START
145         if(TOF){
146             for(tofHit=(AliTOFhit*)TOF->FirstHit(-1);tofHit;
147                 tofHit=(AliTOFhit*)TOF->NextHit()){
148                 hTOF->Fill(tofHit->GetTof());
149             } // end for tofHit
150         } // end if TOF
151
152         if(TRD) {
153             for(trdHit=(AliTRDhit*)TRD->FirstHit(-1);trdHit;
154                 trdHit=(AliTRDhit*)TRD->NextHit()) {
155                 hTRD->Fill((Float_t)(trdHit->GetCharge()));
156             } // end for
157         } // end if TRD
158         if(CASTOR) {
159             for(castorHit=(AliCASTORhit*)CASTOR->FirstHit(-1);castorHit;
160                 castorHit=(AliCASTORhit*)CASTOR->NextHit()) {
161                 hCASTOR->Fill(TMath::Hypot(castorHit->X(),castorHit->Y()));
162             } // end for
163         } // end if CASTOR
164         if(ZDC){
165             for(zdcHit=(AliZDCHit*)ZDC->FirstHit(-1);zdcHit;
166                 zdcHit=(AliZDCHit*)ZDC->NextHit()){
167                 hZDC->Fill(zdcHit->GetEnergy());
168             } // end for zdcdHit
169         } // end if ZDC
170         if(TPC) {        
171           for(tpcHit=(AliTPChit*)TPC->FirstHit(-1);tpcHit;
172               tpcHit=(AliTPChit*)TPC->NextHit()) {
173             hTPC->Fill((Float_t)(tpcHit->fQ));
174           } // end for tpcHit
175         } // end if TPC
176
177     } // end for track
178
179 //Create a canvas, set the view range, show histograms
180     TCanvas *c0 = new TCanvas("c0","Alice Detectors",400,10,600,700);
181     if(opt=="Logy")c0->SetLogy();
182     if(FMD){
183         hFMD->SetFillColor(42);
184         hFMD->Draw();
185         c0->Print("analHitsFMD.ps");
186     } // end if FMD
187     if(ITS){
188         hITS->SetFillColor(42);
189         hITS->Draw();
190         c0->SaveAs("analHitsITS.ps");
191     } // end if ITS
192     if(MUON){
193         hMUON->SetFillColor(42);
194         hMUON->Draw();
195         c0->SaveAs("analHitsMUON.ps");
196     } // end if MUON
197     if(PHOS){
198         hPHOS->SetFillColor(42);
199         hPHOS->Draw();
200         c0->SaveAs("analHitsPHOS.ps");
201     } // end if PHOS
202     if(PMD){
203         hPMD->SetFillColor(42);
204         hPMD->Draw();
205         c0->SaveAs("analHitsPMD.ps");
206     } // end if PMD
207     if(RICH){
208         hRICH->SetFillColor(42);
209         hRICH->Draw();
210         c0->SaveAs("analHitsRICH.ps");
211     } // end if RICH
212     if(START){
213         hSTART->SetFillColor(42);
214         hSTART->Draw();
215         c0->SaveAs("analHitsSTART.ps");
216     } // end if START
217     if(TOF){
218         hTOF->SetFillColor(42);
219         hTOF->Draw();
220         c0->SaveAs("analHitsTOF.ps");
221     } // end if TOF
222     if(TPC){
223         hTPC->SetFillColor(42);
224         hTPC->Draw();
225         c0->SaveAs("analHitsTPC.ps");
226     } // end if TPC
227     if(TRD){
228         hTRD->SetFillColor(42);
229         hTRD->Draw();
230         c0->SaveAs("analHitsTRD.ps");
231     } // end if TRD
232     if(CASTOR){
233         hCASTOR->SetFillColor(42);
234         hCASTOR->Draw();
235         c0->SaveAs("analHitsCASTOR.ps");
236     } // end if TRD
237     if(ZDC){
238         hZDC->SetFillColor(42);
239         hZDC->Draw();
240         c0->SaveAs("analHitsZDC.ps");
241     } // end if ZDC
242
243 // Clean Up
244     /*
245     if(FMD)    delete hFMD;
246     if(ITS)    delete hITS;
247     if(MUON)   delete hMUON;
248     if(PHOS)   delete hPHOS;
249     if(PMD)    delete hPMD;
250     if(RICH)   delete hRICH;
251     if(START)  delete hSTART;
252     if(TOF)    delete hTOF;
253     if(TPC)    delete hTPC;
254     if(TRD)    delete hTRD;
255     if(CASTOR) delete hCASTOR;
256     if(ZDC)    delete hZDC;
257     */
258 }