Minor fixes
[u/mrichter/AliRoot.git] / FMD / scripts / GetMedia.C
CommitLineData
bf000c32 1//____________________________________________________________________
088f8e79 2//
3// $Id$
4//
5// Script that contains a class to get the media where a certain track
6// was created. This is used for background studies.
7//
8// Use the script `Compile.C' to compile this class using ACLic.
9//
10#include <TH2D.h>
11#include <AliFMDHit.h>
12#include <AliFMDInput.h>
13#include <iostream>
14#include <TStyle.h>
15#include <TArrayF.h>
16#include <TArrayI.h>
17#include <AliRun.h>
18#include <TNtuple.h>
55f0ce5b 19#include <AliDetector.h>
20#include <TFile.h>
21#include <TGeoManager.h>
22#include <TGeoVolume.h>
23#include <TGeoMedium.h>
24#include <TParticle.h>
25#include <THStack.h>
9b48326f 26/** @class Media
27 @brief Media description
28 @ingroup FMD_script
29 */
088f8e79 30struct Media : public TNamed
31{
55f0ce5b 32 TArrayI fMeds;
088f8e79 33 TNtuple* fCount;
55f0ce5b 34 TH1D* fHist;
088f8e79 35 Media(const char* name, const char* title)
36 : TNamed(name, title), fMeds(0)
37 {
55f0ce5b 38 fHist = new TH1D(Form("h%s",name), title, 120, -6, 6);
39 // fHist->SetFillStyle(3001);
40 fCount = new TNtuple(GetName(), GetTitle(),
41 "E:DET:RNG:SEC:STR:X:Y:Z:EDEP:PDG");
088f8e79 42 }
55f0ce5b 43 Media(const char* name, TArrayI* a)
44 : TNamed(name, "Media information"), fMeds(a->fN, a->fArray)
088f8e79 45 {
55f0ce5b 46 fHist = new TH1D(Form("h%s",name), GetTitle(), 120, -6, 6);
47 fHist->SetFillStyle(3001);
48 fHist->SetFillColor(fMeds.fN > 0 ? fMeds[0] : 0);
49 fCount = new TNtuple(GetName(), GetTitle(),
50 "E:DET:RNG:SEC:STR:X:Y:Z:EDEP:PDG");
51 std::cout << "Media list " << name << ":" << std::endl;
52 for (Int_t i = 0; i < fMeds.fN; i++) {
53 std::cout << " " << fMeds[i] << std::flush;
54 if (fMeds[i] == 0 && i > 0) break;
088f8e79 55 }
55f0ce5b 56 std::cout << std::endl;
57 }
58 Bool_t HasMedia(Int_t id)
59 {
60 if (fMeds.fN <= 0) return kFALSE;
61 for (Int_t i = 0; i < fMeds.fN; i++) {
62 if (fMeds[i] == id) return kTRUE;
63 if (fMeds[i] == 0 && i > 0) break;
088f8e79 64 }
55f0ce5b 65 return kFALSE;
66 }
67 void Fill(Int_t ev, AliFMDHit* hit)
68 {
69 Float_t x[10];
70 x[0] = ev;
71 x[1] = hit->Detector();
72 x[2] = Int_t(hit->Ring());
73 x[3] = hit->Sector();
74 x[4] = hit->Strip();
75 x[5] = hit->X();
76 x[6] = hit->Y();
77 x[7] = hit->Z();
78 x[8] = hit->Edep();
79 x[9] = hit->Pdg();
80 // return;
81 fCount->Fill(x);
82 // r = TMath::Sqrt(hit->X() * hit->X() + hit->Y() * hit->Y());
83 // if(r!=0) Float_t wgt=1./(2. * 3.1415*r);
84 Double_t r = TMath::Sqrt(hit->X()*hit->X()+hit->Y()*hit->Y());
85 Double_t t = TMath::ATan2(r, hit->Z());
86 Double_t e = -TMath::Log(TMath::Tan(t/2));
87 fHist->Fill(e);
088f8e79 88 }
89};
90
91
bf000c32 92//____________________________________________________________________
9b48326f 93/** @class GetMedia
94 @brief Get media where a particle is produced
95 @code
55f0ce5b 96 aliroot
97 Root> gAlice->Init("$ALICE_ROOT/FMD/Config.C");
98 Root> .x $ALICE_ROOT/FMD/scriptsWriteMedArrays.C
99 Root> gAlice->Run();
100 Root> .q
101 aliroot
9b48326f 102 Root> .L Compile.C
103 Root> Compile("GetMedia.C")
104 Root> GetMedia c
105 Root> c.Run();
106 @endcode
107 @ingroup FMD_script
108 */
088f8e79 109class GetMedia : public AliFMDInputHits
110{
111private:
112 TString fModList;
113 TObjArray fMedia;
114 Media* fOther;
115 Media* fAll;
55f0ce5b 116 Media* fPrim;
088f8e79 117 Int_t fEv;
55f0ce5b 118 THStack* fStack;
088f8e79 119 TFile* fOutput;
120public:
bf000c32 121 //__________________________________________________________________
088f8e79 122 GetMedia(const char* modlist="FMD:ITS:BODY:ABSO:START:PIPE",
123 const char* output="media.root")
55f0ce5b 124 : fModList(modlist)
088f8e79 125 {
55f0ce5b 126 AddLoad(kGeometry);
127 AddLoad(kKinematics);
128
088f8e79 129 fOutput = TFile::Open(output, "RECREATE");
55f0ce5b 130 fStack = new THStack("summed", "Summed hits");
131 fOther = AddMedia("other", "Unknown media");
132 fAll = AddMedia("All", "All media");
133 fPrim = AddMedia("Primary", "Primary particles");
088f8e79 134 fEv = 0;
55f0ce5b 135
136
137 fAll->fHist->SetFillColor(0);
138 fAll->fHist->SetFillStyle(3004);
139 fPrim->fHist->SetFillColor(1);
140 fPrim->fHist->SetFillStyle(3001);
141 fStack->Add(fPrim->fHist);
088f8e79 142 }
bf000c32 143 //__________________________________________________________________
55f0ce5b 144 Media* AddMedia(const char* name, const char* title, TArrayI* a=0)
088f8e79 145 {
088f8e79 146 Media* media = 0;
55f0ce5b 147 if (!a) media = new Media(name, title);
148 else media = new Media(name, a);
149 if (a) {
150 fMedia.Add(media);
151 fStack->Add(media->fHist);
152 }
153 return media;
154 }
155 //__________________________________________________________________
156 Media* FindMedia(Int_t med)
157 {
158 for (Int_t i = 0; i < fMedia.GetEntries(); i++) {
159 Media* media = static_cast<Media*>(fMedia.At(i));
160 if (!media) continue;
161 if (media->HasMedia(med)) return media;
088f8e79 162 }
163 return 0;
164 }
bf000c32 165 //__________________________________________________________________
088f8e79 166 Bool_t Init()
167 {
55f0ce5b 168 Bool_t ret = AliFMDInputHits::Init();
088f8e79 169 if (!gGeoManager) {
55f0ce5b 170 Error("Init", "No geometry defined - make sure you have that");
088f8e79 171 return kFALSE;
172 }
55f0ce5b 173 if (!fRun) {
174 Error("Init", "gAlice not defined");
088f8e79 175 return kFALSE;
176 }
177 Int_t now = 0;
178 Int_t colon;
55f0ce5b 179 TFile* file = TFile::Open("medid.root", "READ");
180 if (!file) {
181 Error("Init", "couldn't open medid.root");
182 return kFALSE;
183 }
184 fOutput->cd();
088f8e79 185 while ((colon = fModList.Index(":", now)) >= 0) {
55f0ce5b 186 TString d(fModList(now, colon-now));
088f8e79 187 now = colon+1;
55f0ce5b 188 TArrayI* a = 0;
189 file->GetObject(d.Data(), a);
190 if (!a) {
191 Warning("Init", "No medium array for %s", d.Data());
192 continue;
193 }
194 AddMedia(d.Data(),0, a);
088f8e79 195 }
55f0ce5b 196 if (now < fModList.Length()) {
197 colon = fModList.Length();
198 TString d(fModList(now, colon-now));
199 TArrayI* a = 0;
200 file->GetObject(d.Data(), a);
201 if (!a)
202 Warning("Init", "No medium array for %s", d.Data());
203 else
204 AddMedia(d.Data(),0, a);
205 }
206 file->Close();
088f8e79 207 if (fMedia.GetEntries() <= 0) return kFALSE;
55f0ce5b 208 return ret;
088f8e79 209 }
bf000c32 210 //__________________________________________________________________
9b48326f 211 /** Begining of event
212 @param ev Event number
213 @return @c false on error */
088f8e79 214 Bool_t Begin(Int_t ev)
215 {
216 fEv = ev;
217 return AliFMDInputHits::Begin(ev);
218 }
bf000c32 219 //__________________________________________________________________
088f8e79 220 Bool_t ProcessHit(AliFMDHit* hit, TParticle* track)
221 {
222 if (!hit || !track) {
223 std::cout << "No hit or track (hit: " << hit << ", track: "
224 << track << ")" << std::endl;
225 return kFALSE;
226 }
55f0ce5b 227 fAll->Fill(fEv, hit);
228 if (track->IsPrimary()) {
229 fPrim->Fill(fEv, hit);
230 return kTRUE;
231 }
232
088f8e79 233 // Get production vertex
234 Double_t vx = track->Vx();
235 Double_t vy = track->Vy();
236 Double_t vz = track->Vz();
088f8e79 237 // Get node
238 TGeoNode* prodNode = gGeoManager->FindNode(vx,vy,vz);
55f0ce5b 239 if (!prodNode) {
240 Warning("ProcessHit", "didn't find a node for production vertex");
241 return kTRUE;
242 }
243 //Info("ProcessHit", "Production vertex in node %s", prodNode->GetName());
244
088f8e79 245 // Get volume
246 TGeoVolume* prodVol = prodNode->GetVolume();
55f0ce5b 247 if (!prodVol) {
248 Warning("ProcessHit", "didn't find a volume for production vertex");
249 return kTRUE;
250 }
251 //Info("ProcessHit", "Production vertex in volume %s",prodVol->GetName());
088f8e79 252
253 // Med medium
254 TGeoMedium* prodMed = prodVol->GetMedium();
55f0ce5b 255 if (!prodMed) {
256 Warning("ProcessHit", "didn't find a medium for production vertex");
257 return kTRUE;
258 }
259 // Info("ProcessHit", "Production vertex in medium %s %d",
260 // prodMed->GetName(), prodMed->GetId());
088f8e79 261
088f8e79 262
55f0ce5b 263 Media* media = FindMedia(prodMed->GetId());
264 if (!media) {
265 Warning("ProcessHit", "Media not found %s (%d)",
266 prodMed->GetName(), prodMed->GetId());
267 media = fOther;
268 }
269 // Info("ProcessHit", "Adding to %s", media->GetName());
270 media->Fill(fEv, hit);
088f8e79 271 return kTRUE;
272 }
bf000c32 273 //__________________________________________________________________
088f8e79 274 Bool_t Finish()
275 {
55f0ce5b 276 fOutput->cd();
277 fStack->Write();
278 fStack->Draw();
088f8e79 279 fOutput->Write();
280 fOutput->Close();
55f0ce5b 281 fOutput = 0;
282 fMedia.Delete();
088f8e79 283 return kTRUE;
284 }
285};
286
287//____________________________________________________________________
288//
289// EOF
290//