Changes for Root6 (Mikolaj)
[u/mrichter/AliRoot.git] / ACORDE / macros / TestACORDEHits.C
1 /////////////////////////////////////////////////////////////////////////
2 //   This macro is a small example of a ROOT macro
3 //   illustrating how to read the hits in ACORDE,
4 //   fill some histograms and store it on a file
5 //
6 //   Author:
7 //
8 //      Eleazar Cuautle <ecuautle@nucleares.unam.mx>
9 //
10 //   Also:
11 //
12 //      Mario Rodriguez Cahuantzi <mrodrigu@mail.cern.ch>
13 //             
14 /////////////////////////////////////////////////////////////////////////
15
16 void TestACORDEHits (const char *filename="galice.root",Int_t numberOfEvents=0){
17    // delete existing gAlice
18   // Dynamically link some shared libs
19   if (gClassTable->GetID("AliRun") < 0) 
20     {
21       gROOT->LoadMacro("loadlibs.C");
22       loadlibs();
23     }
24   if (gAlice)
25     {
26       //delete gAlice->GetRunLoader();
27       delete gAlice;
28       gAlice = 0x0;
29     }
30   
31   AliRunLoader *rl = AliRunLoader::Open("galice.root",
32                                         AliConfig::GetDefaultEventFolderName(),"read");
33   if (!rl)
34     {
35       cerr<<"Can't load RunLoader from file! \n";      return 0x0;
36     }
37   
38   rl->LoadgAlice();
39   gAlice = rl->GetAliRun();
40   if (!gAlice)
41     {
42       cerr << " AliRun object not found on file \n";     return 0x0; 
43     }
44   rl->LoadHeader();
45   
46   // Get the pointer to the ACORDE detector
47   AliLoader *acordel = rl->GetLoader("ACORDELoader");
48   AliACORDE *ACORDE = (AliACORDE *) gAlice->GetDetector("ACORDE");
49   if (ACORDE == 0x0 || acordel == 0x0) {
50     cerr << " Can not find ACORDE or ACORDELoader \n";    return 0x0;  
51   }
52   //Determine how many events would you like to analyze
53   if (numberOfEvents==0) numberOfEvents=(Int_t)(rl->GetNumberOfEvents());
54   
55   //create histograms
56
57   TH1F  *fACORDEEnergy = new TH1F("fACORDEEnergy","Energy distribution of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE); Energy (GeV)",1200,0,360);
58   TH1F  *fACORDEEnergyLoss = new TH1F("fACORDEEnergyLoss","Loss energy distribution  of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE);Energy loss (Gev)",100,0,0.030);
59   TH1F  *fACORDEAzimuth = new TH1F("fACORDEAzimuth","Azimuth angle distribution of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE); Azimuth angle",360,-180,180);
60   TH1F  *fACORDEPolar = new TH1F("fACORDEPolar","Polar angle distribution  of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE); Polar angle",180,0,180);
61   TH2F  *fACORDExz = new TH2F("fACORDExz" ,"Distribution XZ  of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE); X (cms); Z (cms) ",900,-850.,850.,1200,-500.,500.); 
62   TH2F  *fACORDEAzimPol =  new TH2F("fACORDEAzimPol" ,"Azimuth vs Polar angle distribution ",360,-180.,180.,180,0.,180.);
63   TH1F  *fACORDEPx = new TH1F("fACORDEPx" ,"P_{x} distribution  of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE);GeV/c;dN/dP_{x}",60,-120.,120.);
64   TH1F  *fACORDEPy = new TH1F("fACORDEPy" ,"P_{y} distribution  of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE);GeV/c;dN/dP_{y}",60,-120.,120.);
65   TH1F  *fACORDEPz = new TH1F("fACORDEPz" ,"P_{z} distribution  of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE);GeV/c;dN/dP_{z}",60,-120.,120.);
66   TH1F  *fACORDEPt = new TH1F("fACORDEPt" ,"P_{t} distribution  of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE);GeV/c;dN/dP_{t}",1200,0.,360.);
67   TH2F  *fELossAzimuth = new TH2F("fELossAzimuth","Dist. Azimuth Angle VS Energy Loss; Azimuth angle; Energy (GeV)",360,0,350,100,0,0.030); 
68   TH2F  *fEnergyAzimuth = new TH2F("fEnergyAzimuth","Dist. Azimuth Angle VS Energy; Azimuth angle; Energy (GeV)",360,0,360,300,0,350); 
69   TH2F  *fELossPolar = new TH2F("fELossPolar","Dist. Polar Angle VS Energy Loss; Polar angle; Energy (GeV)",180,0,180,100,0,0.030); 
70   TH2F  *fEnergyPolar = new TH2F("fEnergyPolar","Dist. Polar Angle VS Energy; Polar angle; Energy (GeV)",180,0,180,300,0,350); 
71
72  //------------------------------------------------------------------//
73   
74   for (Int_t ievent=0; ievent<numberOfEvents; ievent++) {
75     if ((ievent%10) == 0)  printf ("Processing event %d \n", ievent);
76     rl->GetEvent(ievent);
77     
78     // Get the pointer Hit tree
79     acordel->LoadHits();
80     TTree *hitTree = acordel->TreeH();
81     ACORDE->SetTreeAddress();
82     if (!hitTree) {
83       cout << " No TreeH found" << endl;      return 0x0; //rc;
84     }
85     
86     rl->LoadKinematics();
87     Int_t nTrack = (Int_t) hitTree->GetEntries();
88     //    cout <<"<AliACORDEanalyzeHits> Found "<< nTrack <<" primary particles with hits \n";
89     
90     // Start loop on tracks in the hits containers
91     for(Int_t iTrack=0; iTrack<nTrack;iTrack++){
92       ACORDE->ResetHits();
93       hitTree->GetEvent(iTrack);          
94       if(ACORDE) {       
95         for(acordeHit=(AliACORDEhit*)ACORDE->FirstHit(-1);acordeHit;
96             acordeHit=(AliACORDEhit*)ACORDE->NextHit()) {
97                 fACORDEEnergy->Fill(acordeHit->Energy());
98                 fACORDEEnergyLoss->Fill(acordeHit->Eloss());
99                 fACORDEAzimuth->Fill(acordeHit->AzimuthAngle());
100                 fACORDEPolar->Fill(acordeHit->PolarAngle());
101                 fACORDExz->Fill(acordeHit->X(),acordeHit->Z());
102                 fELossAzimuth->Fill((Float_t)(acordeHit->AzimuthAngle()),(Float_t)(acordeHit->Eloss()));
103                 fACORDEPx->Fill( (Float_t)(acordeHit->Px()));
104                 fACORDEPy->Fill( (Float_t)(acordeHit->Py()));
105                 fACORDEPz->Fill( (Float_t)(acordeHit->Pz()));
106                 fACORDEPt->Fill( TMath::Sqrt( (acordeHit->Px())*(acordeHit->Px())+
107                                              (acordeHit->Py())*(acordeHit->Py())) );
108                 fACORDEAzimPol->Fill( (Float_t)(acordeHit->AzimuthAngle()),
109                                 (Float_t)(acordeHit->PolarAngle()));
110                 fEnergyAzimuth->Fill(acordeHit->Energy(),acordeHit->AzimuthAngle());
111                 fEnergyPolar->Fill(acordeHit->Energy(),acordeHit->PolarAngle());
112                 fELossPolar->Fill(acordeHit->PolarAngle(),acordeHit->Eloss());
113         //cout << "Muon triggered by ACORDE "<<endl;
114         //cout << "Energy: "<<acordeHit->Energy()<<" EnergyLoss: "<<acordeHit->Eloss()<<" Polar angle: "<<acordeHit->PolarAngle()<<" Azimuth angle: "<<acordeHit->AzimuthAngle()<<endl;
115
116         } // end for acordeHits
117       } // end if ACORDE
118     } // end for iTrack  
119   }
120
121  // Put the histograms in a TCanvas
122
123         TCanvas *c1 = new TCanvas("c1","ACORDE Hits distribution of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE)");
124         c1->Divide(2,1);
125
126         c1->cd(1)->SetLogy();
127         fACORDEEnergy->SetMarkerStyle(kFullCircle);
128         fACORDEEnergy->SetMarkerSize(0.7);
129         fACORDEEnergy->Draw("E");
130
131         c1->cd(2)->SetLogy();
132         fACORDEEnergyLoss->SetMarkerStyle(kFullCircle);
133         fACORDEEnergyLoss->SetMarkerSize(0.7);
134         fACORDEEnergyLoss->Draw("E");
135         
136         c1->Draw();
137
138         TCanvas *c2 = new TCanvas("c2","ACORDE Hits distribution of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE)");
139         c2->Divide(2,2);
140         c2->cd(1);
141         fACORDEAzimuth->Draw();
142
143         c2->cd(2);
144         fACORDEPolar->Draw();
145
146         c2->cd(3);
147         gStyle->SetPalette(1);
148         fACORDExz->DrawCopy("colz");
149         
150         c2->cd(4);
151         gStyle->SetPalette(1);
152         fACORDEAzimPol->DrawCopy("colz");
153
154         c2->Draw();
155
156         TCanvas *c3 = new TCanvas("c3","ACORDE Hits distribution of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE)");
157         c3->Divide(2,2);
158
159         c3->cd(1)->SetLogy();
160         fACORDEPx->SetMarkerStyle(kFullCircle);
161         fACORDEPx->SetMarkerSize(0.7);
162         fACORDEPx->Draw("E");
163
164         c3->cd(2)->SetLogy();
165         fACORDEPy->SetMarkerStyle(kFullCircle);
166         fACORDEPy->SetMarkerSize(0.7);
167         fACORDEPy->Draw("E");
168
169         c3->cd(3)->SetLogy();
170         fACORDEPz->SetMarkerStyle(kFullCircle);
171         fACORDEPz->SetMarkerSize(0.7);
172         fACORDEPz->Draw("E");
173
174         c3->cd(4)->SetLogy();
175         fACORDEPt->SetMarkerStyle(kFullCircle);
176         fACORDEPt->SetMarkerSize(0.7);
177         fACORDEPt->Draw("E");
178
179         c3->Draw();
180
181         TCanvas *c4 = new TCanvas("c4","ACORDE Hits distribution of muons generated by AliGenCosmicParam at z = 900 cms (Just above ACORDE)");
182         c4->Divide(2,2);
183
184         c4->cd(1);
185         gStyle->SetPalette(1);
186         fEnergyPolar->DrawCopy("colz");
187
188         c4->cd(2);
189         gStyle->SetPalette(1);
190         fEnergyAzimuth->DrawCopy("colz");
191
192         c4->cd(3);
193         gStyle->SetPalette(1);
194         fELossPolar->DrawCopy("colz");
195
196         c4->cd(4);
197         gStyle->SetPalette(1);
198         fELossAzimuth->DrawCopy("colz");
199         
200         c4->Draw();
201
202
203  // save histos in a root file
204   TFile *fout = new TFile("ACORDE_hits.root","RECREATE");
205         fACORDEEnergy->Write();
206         fACORDEEnergyLoss->Write();
207         fACORDEAzimuth->Write();
208         fACORDEPolar->Write();
209         fACORDExz->Write();
210         fELossAzimuth->Write();
211         fACORDEPx->Write();
212         fACORDEPy->Write();
213         fACORDEPz->Write();
214         fACORDEPt->Write();
215         fACORDEAzimPol->Write();
216         c1->Write();
217         c2->Write();
218         c3->Write();
219         c4->Write();
220 }