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 //