AliRoot Core  3dc7879 (3dc7879)
GetMedia.C
Go to the documentation of this file.
1 //____________________________________________________________________
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>
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>
30 struct Media : public TNamed
31 {
32  TArrayI fMeds;
33  TNtuple* fCount;
34  TH1D* fHist;
35  Media(const char* name, const char* title)
36  : TNamed(name, title), fMeds(0)
37  {
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");
42  }
43  Media(const char* name, TArrayI* a)
44  : TNamed(name, "Media information"), fMeds(a->fN, a->fArray)
45  {
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;
55  }
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;
64  }
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);
88  }
89 };
90 
91 
92 //____________________________________________________________________
109 class GetMedia : public AliFMDInput
110 {
111 private:
112  TString fModList;
117  Int_t fEv;
118  THStack* fStack;
119  TFile* fOutput;
120 public:
121  //__________________________________________________________________
122  GetMedia(const char* modlist="FMD:ITS:BODY:ABSO:T0:PIPE",
123  const char* output="media.root")
124  : fModList(modlist)
125  {
126  AddLoad(kGeometry);
127  AddLoad(kHits);
128  AddLoad(kKinematics);
129 
130  fOutput = TFile::Open(output, "RECREATE");
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");
135  fEv = 0;
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);
143  }
144  //__________________________________________________________________
145  Media* AddMedia(const char* name, const char* title, TArrayI* a=0)
146  {
147  Media* media = 0;
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;
163  }
164  return 0;
165  }
166  //__________________________________________________________________
167  Bool_t Init()
168  {
169  Bool_t ret = AliFMDInputHits::Init();
170  if (!gGeoManager) {
171  Error("Init", "No geometry defined - make sure you have that");
172  return kFALSE;
173  }
174  if (!fRun) {
175  Error("Init", "gAlice not defined");
176  return kFALSE;
177  }
178  Int_t now = 0;
179  Int_t colon;
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();
186  while ((colon = fModList.Index(":", now)) >= 0) {
187  TString d(fModList(now, colon-now));
188  now = colon+1;
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);
196  }
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();
208  if (fMedia.GetEntries() <= 0) return kFALSE;
209  return ret;
210  }
211  //__________________________________________________________________
215  Bool_t Begin(Int_t ev)
216  {
217  fEv = ev;
218  return AliFMDInputHits::Begin(ev);
219  }
220  //__________________________________________________________________
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  }
228  fAll->Fill(fEv, hit);
229  if (track->IsPrimary()) {
230  fPrim->Fill(fEv, hit);
231  return kTRUE;
232  }
233 
234  // Get production vertex
235  Double_t vx = track->Vx();
236  Double_t vy = track->Vy();
237  Double_t vz = track->Vz();
238  // Get node
239  TGeoNode* prodNode = gGeoManager->FindNode(vx,vy,vz);
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 
246  // Get volume
247  TGeoVolume* prodVol = prodNode->GetVolume();
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());
253 
254  // Med medium
255  TGeoMedium* prodMed = prodVol->GetMedium();
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());
262 
263 
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);
272  return kTRUE;
273  }
274  //__________________________________________________________________
275  Bool_t Finish()
276  {
277  fOutput->cd();
278  fStack->Write();
279  fStack->Draw();
280  fOutput->Write();
281  fOutput->Close();
282  fOutput = 0;
283  fMedia.Delete();
284  return kTRUE;
285  }
286 };
287 
288 //____________________________________________________________________
289 //
290 // EOF
291 //
Bool_t Finish()
Definition: GetMedia.C:275
Bool_t Init()
Definition: GetMedia.C:167
Media * fPrim
Definition: GetMedia.C:116
TFile * Open(const char *filename, Long64_t &nevents)
Media * AddMedia(const char *name, const char *title, TArrayI *a=0)
Definition: GetMedia.C:145
Hit in the FMD.
#define TObjArray
Int_t fEv
Definition: GetMedia.C:117
Get media where a particle is produced.
Definition: GetMedia.C:109
FMD utility classes for reading FMD data.
AliRunLoader * Init()
Media(const char *name, TArrayI *a)
Definition: GetMedia.C:43
TFile * fOutput
Definition: GetMedia.C:119
TNtuple * fCount
Definition: GetMedia.C:33
void Fill(Int_t ev, AliFMDHit *hit)
Definition: GetMedia.C:67
AliTPCfastTrack * track
TString fModList
Definition: GetMedia.C:112
Float_t Z() const
Definition: AliHit.h:23
TVectorD vz
Definition: driftITSTPC.C:87
Media(const char *name, const char *title)
Definition: GetMedia.C:35
Float_t X() const
Definition: AliHit.h:21
TH1D * fHist
Definition: GetMedia.C:34
Float_t Y() const
Definition: AliHit.h:22
Bool_t ProcessHit(AliFMDHit *hit, TParticle *track)
Definition: GetMedia.C:221
Float_t Edep() const
Definition: AliFMDHit.h:82
Char_t Ring() const
Definition: AliFMDHit.h:76
Media * fOther
Definition: GetMedia.C:114
Media description.
Definition: GetMedia.C:30
UShort_t Detector() const
Definition: AliFMDHit.h:74
TGeoManager * gGeoManager
Media * FindMedia(Int_t med)
Definition: GetMedia.C:157
Media * fAll
Definition: GetMedia.C:115
UShort_t Sector() const
Definition: AliFMDHit.h:78
Base class for reading in various FMD data. The class loops over all found events. For each event the specified data is read in. The class then loops over all elements of the read data, and process these with user defined code.
Definition: AliFMDInput.h:106
UShort_t Strip() const
Definition: AliFMDHit.h:80
THStack * fStack
Definition: GetMedia.C:118
AliFMDhit is the hit class for the FMD. Hits are the information that comes from a Monte Carlo at eac...
Definition: AliFMDHit.h:30
GetMedia(const char *modlist="FMD:ITS:BODY:ABSO:T0:PIPE", const char *output="media.root")
Definition: GetMedia.C:122
Int_t Pdg() const
Definition: AliFMDHit.h:96
TObjArray fMedia
Definition: GetMedia.C:113
Bool_t HasMedia(Int_t id)
Definition: GetMedia.C:58
TArrayI fMeds
Definition: GetMedia.C:32
Bool_t Begin(Int_t ev)
Definition: GetMedia.C:215