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 AliFMDInput
 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: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   //__________________________________________________________________
 212   /** Begining of event
 213       @param ev Event number
 214       @return @c false on error */
 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 //