Fixed some coding style violations.
[u/mrichter/AliRoot.git] / FMD / scripts / GetMedia.C
1 //
2 // $Id$
3 //
4 // Script that contains a class to get the media where a certain track
5 // was created.   This is used for background studies. 
6 //
7 // Use the script `Compile.C' to compile this class using ACLic. 
8 //
9 #include <TH2D.h>
10 #include <AliFMDHit.h>
11 #include <AliFMDInput.h>
12 #include <iostream>
13 #include <TStyle.h>
14 #include <TArrayF.h>
15 #include <TArrayI.h>
16 #include <AliRun.h>
17 #include <TNtuple.h>
18
19 struct Media : public TNamed
20 {
21   TArrayI*        fMeds;
22   TNtuple*        fCount;
23   Media(const char* name, const char* title) 
24     : TNamed(name, title), fMeds(0) 
25   {
26     fCount = new TNtuple(GetName(), "E/I:DET:RNG:SEC:STR:X/F:Y:Z:EDEP:PDG/I");
27   }
28   Media(const char* name) 
29     : TNamed(name, "Media information"), fMeds(0), fCount(0)
30   {
31     AliDetector* det = gAlice->GetDetector(name);
32     if (!det) {
33       Warning("Media", "Detector %s not found in gAlice", name);
34       return;
35     }
36     fMeds = det->GetIdtmed();
37     if (!fMeds) {
38       Warning("Media", "No mediums for detector %s found", name);
39       return;
40     }
41     fCount = new TNtuple(GetName(), "E/I:DET:RNG:SEC:STR:X/F:Y:Z:EDEP:PDG/I");
42   }
43 };
44
45
46 class GetMedia : public AliFMDInputHits
47 {
48 private:
49   TString    fModList;
50   TObjArray  fMedia;
51   Media*     fOther;
52   Media*     fAll;
53   Int_t      fEv;
54   TFile*     fOutput;
55 public:
56   GetMedia(const char* modlist="FMD:ITS:BODY:ABSO:START:PIPE", 
57            const char* output="media.root") 
58   { 
59     fOutput = TFile::Open(output, "RECREATE");
60     fOther  = new Media("other", "Unknown media"),
61     fAll    = new Media("All", "All media")
62     fEv     = 0;
63     AddLoad(kKinematics);
64   }
65   Media* FindMedia(Int_t med) 
66   {
67     TIter next(&fMedia);
68     Media* media = 0;
69     while ((media == static_cast<Media*>(next()))) {
70       if (!media->fMeds) continue;
71       TArrayI& ids = *(media->fMeds);
72       for (Int_t i = 0; i < ids.fN; i++) 
73         if (med == ids[i]) return media;
74     }
75     return 0;
76   }
77   Bool_t Init()
78   {
79     if (!gGeoManager) {
80       Error("GetMedia", "No geometry defined - make sure you have that");
81       return kFALSE;
82     }
83     if (!gAlice) {
84       Error("GetMedia", "gAlice not defined");
85       return kFALSE;
86     }
87     Int_t now = 0;
88     Int_t colon;
89     while ((colon = fModList.Index(":", now)) >= 0) {
90       fMedia.Add(new Media(fModList(now, colon-now)));
91       now = colon+1;
92     }
93     if (now < fModList.Length()) 
94       fMedia.Add(new Media(now, fModList.Length()-now));
95     if (fMedia.GetEntries() <= 0) return kFALSE;
96     return AliFMDInputHits::Init();
97   }  
98   Bool_t Begin(Int_t ev) 
99   {
100     fEv = ev;
101     return AliFMDInputHits::Begin(ev);
102   }
103   Bool_t ProcessHit(AliFMDHit* hit, TParticle* track) 
104   {
105     if (!hit || !track) {
106       std::cout << "No hit or track (hit: " << hit << ", track: " 
107                 << track << ")" << std::endl;
108       return kFALSE;
109     }
110     // Get production vertex 
111     Double_t vx = track->Vx();
112     Double_t vy = track->Vy();
113     Double_t vz = track->Vz();
114     
115     fAll->fCounts->Fill(fEv, hit->Detector(),Int_t(hit->Ring()),
116                         hit->Sector(), hit->Strip(), 
117                         hit->X(), hit->Y(), hit->Z(), hit->Edep(), 
118                         hit->Pdg());
119     // Get node 
120     TGeoNode* prodNode = gGeoManager->FindNode(vx,vy,vz);
121     if (!prodNode) return kTRUE;
122     
123     // Get volume 
124     TGeoVolume* prodVol = prodNode->GetVolume();
125     if (!prodVol) return kTRUE;
126     
127     // Med medium 
128     TGeoMedium* prodMed = prodVol->GetMedium();
129     if (!prodMed) return kTRUE;
130     
131     Media* media = FindMedia(prodMed->GetUniqueID());
132     if (media) media = fOther; 
133     
134     // r = TMath::Sqrt(hit->X() * hit->X() + hit->Y() * hit->Y());
135     // if(r!=0) Float_t wgt=1./(2. * 3.1415*r);
136     media->fCounts->Fill(fEv, hit->Detector(),Int_t(hit->Ring()),
137                          hit->Sector(), hit->Strip(), 
138                          hit->X(), hit->Y(), hit->Z(), hit->Edep(), 
139                          hit->Pdg());
140     return kTRUE;
141   }
142   Bool_t Finish()
143   {
144     fOutput->Write();
145     fOutput->Close();
146     return kTRUE;
147   }
148 };
149
150 //____________________________________________________________________
151 //
152 // EOF
153 //