A lot of changes after detector review:
[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 */
8ec606c2 109class GetMedia : public AliFMDInput
088f8e79 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 //__________________________________________________________________
ababa197 122 GetMedia(const char* modlist="FMD:ITS:BODY:ABSO:T0:PIPE",
088f8e79 123 const char* output="media.root")
55f0ce5b 124 : fModList(modlist)
088f8e79 125 {
55f0ce5b 126 AddLoad(kGeometry);
8ec606c2 127 AddLoad(kHits);
55f0ce5b 128 AddLoad(kKinematics);
129
088f8e79 130 fOutput = TFile::Open(output, "RECREATE");
55f0ce5b 131 fStack = new THStack("summed", "Summed hits");
132 fOther = AddMedia("other", "Unknown media");
133 fAll = AddMedia("All", "All media");
134 fPrim = AddMedia("Primary", "Primary particles");
088f8e79 135 fEv = 0;
55f0ce5b 136
137
138 fAll->fHist->SetFillColor(0);
139 fAll->fHist->SetFillStyle(3004);
140 fPrim->fHist->SetFillColor(1);
141 fPrim->fHist->SetFillStyle(3001);
142 fStack->Add(fPrim->fHist);
088f8e79 143 }
bf000c32 144 //__________________________________________________________________
55f0ce5b 145 Media* AddMedia(const char* name, const char* title, TArrayI* a=0)
088f8e79 146 {
088f8e79 147 Media* media = 0;
55f0ce5b 148 if (!a) media = new Media(name, title);
149 else media = new Media(name, a);
150 if (a) {
151 fMedia.Add(media);
152 fStack->Add(media->fHist);
153 }
154 return media;
155 }
156 //__________________________________________________________________
157 Media* FindMedia(Int_t med)
158 {
159 for (Int_t i = 0; i < fMedia.GetEntries(); i++) {
160 Media* media = static_cast<Media*>(fMedia.At(i));
161 if (!media) continue;
162 if (media->HasMedia(med)) return media;
088f8e79 163 }
164 return 0;
165 }
bf000c32 166 //__________________________________________________________________
088f8e79 167 Bool_t Init()
168 {
55f0ce5b 169 Bool_t ret = AliFMDInputHits::Init();
088f8e79 170 if (!gGeoManager) {
55f0ce5b 171 Error("Init", "No geometry defined - make sure you have that");
088f8e79 172 return kFALSE;
173 }
55f0ce5b 174 if (!fRun) {
175 Error("Init", "gAlice not defined");
088f8e79 176 return kFALSE;
177 }
178 Int_t now = 0;
179 Int_t colon;
55f0ce5b 180 TFile* file = TFile::Open("medid.root", "READ");
181 if (!file) {
182 Error("Init", "couldn't open medid.root");
183 return kFALSE;
184 }
185 fOutput->cd();
088f8e79 186 while ((colon = fModList.Index(":", now)) >= 0) {
55f0ce5b 187 TString d(fModList(now, colon-now));
088f8e79 188 now = colon+1;
55f0ce5b 189 TArrayI* a = 0;
190 file->GetObject(d.Data(), a);
191 if (!a) {
192 Warning("Init", "No medium array for %s", d.Data());
193 continue;
194 }
195 AddMedia(d.Data(),0, a);
088f8e79 196 }
55f0ce5b 197 if (now < fModList.Length()) {
198 colon = fModList.Length();
199 TString d(fModList(now, colon-now));
200 TArrayI* a = 0;
201 file->GetObject(d.Data(), a);
202 if (!a)
203 Warning("Init", "No medium array for %s", d.Data());
204 else
205 AddMedia(d.Data(),0, a);
206 }
207 file->Close();
088f8e79 208 if (fMedia.GetEntries() <= 0) return kFALSE;
55f0ce5b 209 return ret;
088f8e79 210 }
bf000c32 211 //__________________________________________________________________
9b48326f 212 /** Begining of event
213 @param ev Event number
214 @return @c false on error */
088f8e79 215 Bool_t Begin(Int_t ev)
216 {
217 fEv = ev;
218 return AliFMDInputHits::Begin(ev);
219 }
bf000c32 220 //__________________________________________________________________
088f8e79 221 Bool_t ProcessHit(AliFMDHit* hit, TParticle* track)
222 {
223 if (!hit || !track) {
224 std::cout << "No hit or track (hit: " << hit << ", track: "
225 << track << ")" << std::endl;
226 return kFALSE;
227 }
55f0ce5b 228 fAll->Fill(fEv, hit);
229 if (track->IsPrimary()) {
230 fPrim->Fill(fEv, hit);
231 return kTRUE;
232 }
233
088f8e79 234 // Get production vertex
235 Double_t vx = track->Vx();
236 Double_t vy = track->Vy();
237 Double_t vz = track->Vz();
088f8e79 238 // Get node
239 TGeoNode* prodNode = gGeoManager->FindNode(vx,vy,vz);
55f0ce5b 240 if (!prodNode) {
241 Warning("ProcessHit", "didn't find a node for production vertex");
242 return kTRUE;
243 }
244 //Info("ProcessHit", "Production vertex in node %s", prodNode->GetName());
245
088f8e79 246 // Get volume
247 TGeoVolume* prodVol = prodNode->GetVolume();
55f0ce5b 248 if (!prodVol) {
249 Warning("ProcessHit", "didn't find a volume for production vertex");
250 return kTRUE;
251 }
252 //Info("ProcessHit", "Production vertex in volume %s",prodVol->GetName());
088f8e79 253
254 // Med medium
255 TGeoMedium* prodMed = prodVol->GetMedium();
55f0ce5b 256 if (!prodMed) {
257 Warning("ProcessHit", "didn't find a medium for production vertex");
258 return kTRUE;
259 }
260 // Info("ProcessHit", "Production vertex in medium %s %d",
261 // prodMed->GetName(), prodMed->GetId());
088f8e79 262
088f8e79 263
55f0ce5b 264 Media* media = FindMedia(prodMed->GetId());
265 if (!media) {
266 Warning("ProcessHit", "Media not found %s (%d)",
267 prodMed->GetName(), prodMed->GetId());
268 media = fOther;
269 }
270 // Info("ProcessHit", "Adding to %s", media->GetName());
271 media->Fill(fEv, hit);
088f8e79 272 return kTRUE;
273 }
bf000c32 274 //__________________________________________________________________
088f8e79 275 Bool_t Finish()
276 {
55f0ce5b 277 fOutput->cd();
278 fStack->Write();
279 fStack->Draw();
088f8e79 280 fOutput->Write();
281 fOutput->Close();
55f0ce5b 282 fOutput = 0;
283 fMedia.Delete();
088f8e79 284 return kTRUE;
285 }
286};
287
288//____________________________________________________________________
289//
290// EOF
291//