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
  20 /** @class Media
  21     @brief Media description
  22     @ingroup FMD_script
  23  */
  24 struct Media : public TNamed
  25 {
  26   TArrayI*        fMeds;
  27   TNtuple*        fCount;
  28   Media(const char* name, const char* title)
  29     : TNamed(name, title), fMeds(0)
  30   {
  31     fCount = new TNtuple(GetName(), "E/I:DET:RNG:SEC:STR:X/F:Y:Z:EDEP:PDG/I");
  32   }
  33   Media(const char* name)
  34     : TNamed(name, "Media information"), fMeds(0), fCount(0)
  35   {
  36     AliDetector* det = gAlice->GetDetector(name);
  37     if (!det) {
  38       Warning("Media", "Detector %s not found in gAlice", name);
  39       return;
  40     }
  41     fMeds = det->GetIdtmed();
  42     if (!fMeds) {
  43       Warning("Media", "No mediums for detector %s found", name);
  44       return;
  45     }
  46     fCount = new TNtuple(GetName(), "E/I:DET:RNG:SEC:STR:X/F:Y:Z:EDEP:PDG/I");
  47   }
  48 };
  49
  50
  51 //____________________________________________________________________
  52 /** @class GetMedia
  53     @brief Get media where a particle is produced
  54     @code
  55     Root> .L Compile.C
  56     Root> Compile("GetMedia.C")
  57     Root> GetMedia c
  58     Root> c.Run();
  59     @endcode
  60     @ingroup FMD_script
  61  */
  62 class GetMedia : public AliFMDInputHits
  63 {
  64 private:
  65   TString    fModList;
  66   TObjArray  fMedia;
  67   Media*     fOther;
  68   Media*     fAll;
  69   Int_t      fEv;
  70   TFile*     fOutput;
  71 public:
  72   //__________________________________________________________________
  73   GetMedia(const char* modlist="FMD:ITS:BODY:ABSO:START:PIPE",
  74            const char* output="media.root")
  75   {
  76     fOutput = TFile::Open(output, "RECREATE");
  77     fOther  = new Media("other", "Unknown media"),
  78     fAll    = new Media("All", "All media")
  79     fEv     = 0;
  80     AddLoad(kKinematics);
  81   }
  82   //__________________________________________________________________
  83   Media* FindMedia(Int_t med)
  84   {
  85     TIter next(&fMedia);
  86     Media* media = 0;
  87     while ((media == static_cast<Media*>(next()))) {
  88       if (!media->fMeds) continue;
  89       TArrayI& ids = *(media->fMeds);
  90       for (Int_t i = 0; i < ids.fN; i++)
  91         if (med == ids[i]) return media;
  92     }
  93     return 0;
  94   }
  95   //__________________________________________________________________
  96   Bool_t Init()
  97   {
  98     if (!gGeoManager) {
  99       Error("GetMedia", "No geometry defined - make sure you have that");
 100       return kFALSE;
 101     }
 102     if (!gAlice) {
 103       Error("GetMedia", "gAlice not defined");
 104       return kFALSE;
 105     }
 106     Int_t now = 0;
 107     Int_t colon;
 108     while ((colon = fModList.Index(":", now)) >= 0) {
 109       fMedia.Add(new Media(fModList(now, colon-now)));
 110       now = colon+1;
 111     }
 112     if (now < fModList.Length())
 113       fMedia.Add(new Media(now, fModList.Length()-now));
 114     if (fMedia.GetEntries() <= 0) return kFALSE;
 115     return AliFMDInputHits::Init();
 116   }
 117   //__________________________________________________________________
 118   /** Begining of event
 119       @param ev Event number
 120       @return @c false on error */
 121   Bool_t Begin(Int_t ev)
 122   {
 123     fEv = ev;
 124     return AliFMDInputHits::Begin(ev);
 125   }
 126   //__________________________________________________________________
 127   Bool_t ProcessHit(AliFMDHit* hit, TParticle* track)
 128   {
 129     if (!hit || !track) {
 130       std::cout << "No hit or track (hit: " << hit << ", track: "
 131                 << track << ")" << std::endl;
 132       return kFALSE;
 133     }
 134     // Get production vertex
 135     Double_t vx = track->Vx();
 136     Double_t vy = track->Vy();
 137     Double_t vz = track->Vz();
 138
 139     fAll->fCounts->Fill(fEv, hit->Detector(),Int_t(hit->Ring()),
 140                         hit->Sector(), hit->Strip(),
 141                         hit->X(), hit->Y(), hit->Z(), hit->Edep(),
 142                         hit->Pdg());
 143     // Get node
 144     TGeoNode* prodNode = gGeoManager->FindNode(vx,vy,vz);
 145     if (!prodNode) return kTRUE;
 146
 147     // Get volume
 148     TGeoVolume* prodVol = prodNode->GetVolume();
 149     if (!prodVol) return kTRUE;
 150
 151     // Med medium
 152     TGeoMedium* prodMed = prodVol->GetMedium();
 153     if (!prodMed) return kTRUE;
 154
 155     Media* media = FindMedia(prodMed->GetUniqueID());
 156     if (media) media = fOther;
 157
 158     // r = TMath::Sqrt(hit->X() * hit->X() + hit->Y() * hit->Y());
 159     // if(r!=0) Float_t wgt=1./(2. * 3.1415*r);
 160     media->fCounts->Fill(fEv, hit->Detector(),Int_t(hit->Ring()),
 161                          hit->Sector(), hit->Strip(),
 162                          hit->X(), hit->Y(), hit->Z(), hit->Edep(),
 163                          hit->Pdg());
 164     return kTRUE;
 165   }
 166   //__________________________________________________________________
 167   Bool_t Finish()
 168   {
 169     fOutput->Write();
 170     fOutput->Close();
 171     return kTRUE;
 172   }
 173 };
 174
 175 //____________________________________________________________________
 176 //
 177 // EOF
 178 //