Minor fixes
[u/mrichter/AliRoot.git] / FMD / scripts / GetMedia.C
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>
26 /** @class Media 
27     @brief Media description 
28     @ingroup FMD_script
29  */
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 //____________________________________________________________________
93 /** @class GetMedia
94     @brief Get media where a particle is produced
95     @code 
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
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  */
109 class GetMedia : public AliFMDInputHits
110 {
111 private:
112   TString    fModList;
113   TObjArray  fMedia;
114   Media*     fOther;
115   Media*     fAll;
116   Media*     fPrim;
117   Int_t      fEv;
118   THStack*   fStack;
119   TFile*     fOutput;
120 public:
121   //__________________________________________________________________
122   GetMedia(const char* modlist="FMD:ITS:BODY:ABSO:START:PIPE", 
123            const char* output="media.root") 
124     :  fModList(modlist)
125   { 
126     AddLoad(kGeometry);
127     AddLoad(kKinematics);
128
129     fOutput = TFile::Open(output, "RECREATE");
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");
134     fEv     = 0;
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);
142   }
143   //__________________________________________________________________
144   Media* AddMedia(const char* name, const char* title, TArrayI* a=0)
145   {
146     Media* media = 0;
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;
162     }
163     return 0;
164   }
165   //__________________________________________________________________
166   Bool_t Init()
167   {
168     Bool_t ret = AliFMDInputHits::Init();
169     if (!gGeoManager) {
170       Error("Init", "No geometry defined - make sure you have that");
171       return kFALSE;
172     }
173     if (!fRun) {
174       Error("Init", "gAlice not defined");
175       return kFALSE;
176     }
177     Int_t now = 0;
178     Int_t colon;
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();
185     while ((colon = fModList.Index(":", now)) >= 0) {
186       TString d(fModList(now, colon-now));
187       now = colon+1;
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);
195     }
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();
207     if (fMedia.GetEntries() <= 0) return kFALSE;
208     return ret;
209   }  
210   //__________________________________________________________________
211   /** Begining of event
212       @param ev Event number
213       @return @c false on error */
214   Bool_t Begin(Int_t ev) 
215   {
216     fEv = ev;
217     return AliFMDInputHits::Begin(ev);
218   }
219   //__________________________________________________________________
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     }
227     fAll->Fill(fEv, hit);
228     if (track->IsPrimary()) { 
229       fPrim->Fill(fEv, hit);
230       return kTRUE;
231     }
232       
233     // Get production vertex 
234     Double_t vx = track->Vx();
235     Double_t vy = track->Vy();
236     Double_t vz = track->Vz();
237     // Get node 
238     TGeoNode* prodNode = gGeoManager->FindNode(vx,vy,vz);
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
245     // Get volume 
246     TGeoVolume* prodVol = prodNode->GetVolume();
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());
252     
253     // Med medium 
254     TGeoMedium* prodMed = prodVol->GetMedium();
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());
261     
262     
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);
271     return kTRUE;
272   }
273   //__________________________________________________________________
274   Bool_t Finish()
275   {
276     fOutput->cd();
277     fStack->Write();
278     fStack->Draw();
279     fOutput->Write();
280     fOutput->Close();
281     fOutput = 0;
282     fMedia.Delete();
283     return kTRUE;
284   }
285 };
286
287 //____________________________________________________________________
288 //
289 // EOF
290 //