AliGeant4/macro/hits.C

   1 // $Id$
   2
   3 #include <iostream.h>
   4
   5 void hits()
   6 {
   7   hits("FMD","1");
   8   hits("PHOS","1");
   9   hits("START","1");
  10   hits("TOF","2");
  11   hits("TPC","2");
  12
  13   //hits("ITS","5");
  14   //hits("MUON","0");
  15   //hits("PMD","1");
  16   //hits("RICH","1");
  17 }
  18
  19 void hits(const TString& detName, const TString& detVersion)
  20 {
  21   cout << detName << endl;
  22
  23   // labels
  24   TString g3 = "G3: ";
  25   TString g4 = "G4: ";
  26
  27   // construct file names
  28   TString top = getenv("ALICE_ROOT");
  29   TString f1NameEnd = "test10.root";
  30   TString f2NameEnd = "test20.root";
  31
  32   TString g3File1Name = top + "/test/" + f1NameEnd;
  33   TString g3File2Name = top + "/test/" + f2NameEnd;
  34   TString g4File1Name = top + "/AliGeant4/test/" + f1NameEnd;
  35   TString g4File2Name = top + "/AliGeant4/test/" + f2NameEnd;
  36
  37   cout << "Files: " << endl;
  38   cout << g3File1Name << endl;
  39   cout << g3File2Name << endl;
  40   cout << g4File1Name << endl;
  41   cout << g4File2Name << endl;
  42
  43   // link shared libs
  44   if (gClassTable->GetID("AliRun") < 0) {
  45     gROOT->LoadMacro("loadlibs.C");
  46     loadlibs();
  47   }
  48   cout << "AliRoot libraries were loaded." << endl;
  49
  50   // set histogram ranges
  51   Int_t nbin;
  52   Int_t g3xmin; Int_t g4xmin; Int_t g3xmax; Int_t g4xmax;
  53   Int_t g3ymin; Int_t g4ymin; Int_t g3ymax; Int_t g4ymax;
  54   Int_t g3zmin; Int_t g4zmin; Int_t g3zmax; Int_t g4zmax;
  55   SetHistogramRanges(detName, nbin,
  56                      g3xmin, g3xmax, g4xmin, g4xmax,
  57                      g3ymin, g3ymax, g4ymin, g4ymax,
  58                      g3zmin, g3zmax, g4zmin, g4zmax);
  59
  60
  61   // create histograms
  62   TH1F* g3x  = new TH1F("g3x",  g3 + detName + " hits per x", nbin, g3xmin, g3xmax);
  63   TH1F* g3xh = new TH1F("g3xh", g3 + detName + " hits per x", nbin, g3xmin, g3xmax);
  64   TH1F* g4x  = new TH1F("g4x",  g4 + detName + " hits per x", nbin, g4xmin, g4xmax);
  65   TH1F* g4xh = new TH1F("g4xh", g4 + detName + " hits per x", nbin, g4xmin, g4xmax);
  66
  67   TH1F* g3z  = new TH1F("g3z",  g3 + detName + " hits per z", nbin, g3zmin, g3zmax);
  68   TH1F* g3zh = new TH1F("g3zh", g3 + detName + " hits per z", nbin, g3zmin, g3zmax);
  69   TH1F* g4z  = new TH1F("g4z",  g4 + detName + " hits per z", nbin, g4zmin, g4zmax);
  70   TH1F* g4zh = new TH1F("g4zh", g4 + detName + " hits per z", nbin, g4zmin, g4zmax);
  71
  72   TH1F* g3y  = 0;
  73   TH1F* g3yh = 0;
  74   TH1F* g4y  = 0;
  75   TH1F* g4yh = 0;
  76   if (detName == "PHOS") {
  77     TH1F* g3y  = new TH1F("g3y",  g3 + detName + " hits per y", nbin, g3ymin, g3ymax);
  78     TH1F* g3yh = new TH1F("g3yh", g3 + detName + " hits per y", nbin, g3ymin, g3ymax);
  79     TH1F* g4y  = new TH1F("g4y",  g4 + detName + " hits per y", nbin, g4ymin, g4ymax);
  80     TH1F* g4yh = new TH1F("g4yh", g4 + detName + " hits per y", nbin, g4ymin, g4ymax);
  81   }
  82
  83   cout << "Histograms were created." << endl;
  84
  85   // fill histograms
  86   AliDetector* detector;
  87
  88   detector = LoadDetector(g3File1Name, detName, g3);
  89   FillHistogram(detector, g3, g3x, g3y, g3z);
  90
  91   detector = LoadDetector(g3File2Name, detName,  g3);
  92   FillHistogram(detector, g3, g3xh, g3yh, g3zh);
  93   //TH1F* g3xh = LoadHistograms("g3xh");
  94   //TH1F* g3zh = LoadHistograms("g3zh");
  95
  96   detector = LoadDetector(g4File1Name, detName, g4);
  97   FillHistogram(detector, g4, g4x, g4y, g4z);
  98
  99   detector = LoadDetector(g4File2Name, detName, g4);
 100   FillHistogram(detector, g4, g4xh, g4yh, g4zh);
 101
 102   // compose picture name
 103   TString gifNameBase =  "hits" + detName + "v" + detVersion;
 104   TString title = detName + " hits";
 105
 106   // draw histograms (G3, G4 in one canvas)
 107   DrawHistograms(title, gifNameBase + ".gif", g3xh, g4xh, g3zh, g4zh);
 108   if (detName == "PHOS")
 109     DrawHistograms(title, gifNameBase + "2.gif", g3yh, g4yh, g3zh, g4zh);
 110
 111   //DrawHistograms(title, gifNameBase + "_x.gif", g3x, g4x, g3xh, g4xh);
 112   //DrawHistograms(title, gifNameBase + "_y.gif", g3y, g4y, g3yh, g4yh);
 113   //DrawHistograms(title, gifNameBase + "_z.gif", g3z, g4z, g3zh, g4zh);
 114
 115   // draw histograms (one by one)
 116   DrawHistogram(title, "g3x_"+ gifNameBase + ".gif", g3xh);
 117   DrawHistogram(title, "g4x_"+ gifNameBase + ".gif", g4xh);
 118   DrawHistogram(title, "g3z_"+ gifNameBase + ".gif", g3zh);
 119   DrawHistogram(title, "g4z_"+ gifNameBase + ".gif", g4zh);
 120
 121 }
 122
 123 void SetHistogramRanges(TString& detName, Int_t& nbin,
 124                         Int_t& g3xmin, Int_t& g3xmax,
 125                         Int_t& g4xmin, Int_t& g4xmax,
 126                         Int_t& g3ymin, Int_t& g3ymax,
 127                         Int_t& g4ymin, Int_t& g4ymax,
 128                         Int_t& g3zmin, Int_t& g3zmax,
 129                         Int_t& g4zmin, Int_t& g4zmax)
 130 {
 131   nbin = 200;
 132
 133   g3xmin = 0; g4xmin = g3xmin;
 134   g3xmax = 1; g4xmax = g3xmax;
 135   g3ymin = 0; g4ymin = g3xmin;
 136   g3ymax = 1; g4ymax = g3xmax;
 137   g3zmin = 0; g4zmin = g3zmin;
 138   g3zmax = 1; g4zmax = g3zmax;
 139
 140   if (detName == "FMD") {
 141     g3xmin = -50; g4xmin = g3xmin;
 142     g3xmax =  50; g4xmax = g3xmax;
 143     g3zmin = -700; g4zmin = g3zmin;
 144     g3zmax =  150; g4zmax = g3zmax;
 145   }
 146   else if (detName == "ITS") {
 147     g3xmin = -60; g4xmin = g3xmin;
 148     g3xmax =  60; g4xmax = g3xmax;
 149     g3zmin = -60; g4zmin = g3zmin;
 150     g3zmax =  60; g4zmax = g3zmax;
 151   }
 152   else if (detName == "MUON") {
 153     g3xmin = -450; g4xmin = g3xmin;
 154     g3xmax =  450; g4xmax = g3xmax;
 155     g3zmin =  400; g4zmin = g3zmin;
 156     g3zmax = 1800; g4zmax = g3zmax;
 157   }
 158   else if (detName == "PHOS") {
 159     g3xmin = -400; g4xmin = g3xmin;
 160     g3xmax =  400; g4xmax = g3xmax;
 161     g3ymin =  250; g4ymin = g3ymin;
 162     g3ymax =  500; g4ymax = g3ymax;
 163     g3zmin =  -80; g4zmin = g3zmin;
 164     g3zmax =   80; g4zmax = g3zmax;
 165   }
 166   else if (detName == "PMD") {
 167     g3xmin = -200; g4xmin = g3xmin*10;
 168     g3xmax =  200; g4xmax = g3xmax*10;
 169     g3zmin = -582; g4zmin = g3zmin*10;
 170     g3zmax = -578; g4zmax = g3zmax*10;
 171   }
 172   else if (detName == "RICH") {
 173     g3xmin = -250; g4xmin = g3xmin;
 174     g3xmax =  250; g4xmax = g3xmax;
 175     g3zmin = -250; g4zmin = g3zmin;
 176     g3zmax =  250; g4zmax = g3zmax;
 177   }
 178   else if (detName == "START") {
 179     g3xmin = -10; g4xmin = g3xmin;
 180     g3xmax =  10; g4xmax = g3xmax;
 181     g3zmin = -80; g4zmin = g3zmin;
 182     g3zmax =  80; g4zmax = g3zmax;
 183   }
 184   else if (detName == "TOF") {
 185     g3xmin = -400; g4xmin = g3xmin;
 186     g3xmax =  400; g4xmax = g3xmax;
 187     g3zmin = -400; g4zmin = g3zmin;
 188     g3zmax =  400; g4zmax = g3zmax;
 189   }
 190   else if (detName == "TPC") {
 191     g3xmin = -300; g4xmin = g3xmin;
 192     g3xmax =  300; g4xmax = g3xmax;
 193     g3zmin = -300; g4zmin = g3zmin;
 194     g3zmax =  300; g4zmax = g3zmax;
 195   }
 196 }
 197
 198 AliDetector* LoadDetector(TString& fileName, TString& detName, TString& label)
 199 {
 200   // connect the Root files
 201   TFile* file = new TFile(fileName);
 202   if (!file) {
 203     cerr << "Root file was not found." << endl;
 204     return 0;
 205   }
 206
 207   // get AliRun object from file
 208   if (gAlice) delete gAlice;
 209   gAlice = 0;
 210   gAlice = (AliRun*)file->Get("gAlice");
 211   if (!gAlice) {
 212     cerr << label << "AliRun object not found in a file." << endl;
 213     return 0;
 214   }
 215
 216   // import the hits trees
 217   Int_t nparticles = gAlice->GetEvent(0);
 218   if (nparticles <= 0) {
 219     cerr << label << "No particles in a file." << endl;
 220     return 0;
 221   }
 222   cout << label << "got nparticles = " << nparticles << endl;
 223
 224   // get detector
 225   AliDetector* detector = gAlice->GetDetector(detName);
 226   if (!detector) {
 227     cerr << label << "No detector " << detName << " in a file." << endl;
 228     return 0;
 229   }
 230
 231   return detector;
 232 }
 233
 234 void FillHistogram(AliDetector* detector, TString& label,
 235                    TH1F* hx, TH1F* hy, TH1F* hz)
 236 {
 237   Int_t nofHits = 0;
 238   // get number of primary tracks
 239   Int_t ntracks = gAlice->TreeH()->GetEntries();
 240   cout << label << "got ntracks = " << ntracks << endl;
 241
 242   // loop on tracks in the hits container
 243   for (Int_t i=0; i<ntracks; i++) {
 244     // loop on hits
 245     for(AliHit* hit= detector->FirstHit(i); hit; hit=detector->NextHit()) {
 246
 247       TString detName = detector->GetName();
 248       if (detName == "PHOS") {
 249         // PHOS special
 250         FillPHOSHit(detector, hit, hx, hy, hz);
 251       }
 252       else {
 253         FillHit(hit, hx, hy, hz);
 254       }
 255
 256       nofHits++;
 257     }
 258   }
 259 /*
 260   TFile* file = new TFile("tmp.root","recreate");
 261   hx->Write();
 262   hz->Write();
 263   file->Close();
 264 */
 265   cout << label << "filled " << nofHits << " hits" << endl;
 266 }
 267
 268 TH1F* LoadHistograms(TString name)
 269 {
 270   TFile* file = new TFile("tmp.root");
 271   return (TH1F*)file->Get(name);
 272 }
 273
 274 void FillHit(AliHit* hit, TH1F* hx, TH1F* hy, TH1F* hz)
 275 {
 276   if (hx) hx->Fill(hit->X());
 277   if (hy) hy->Fill(hit->Y());
 278   if (hz) hz->Fill(hit->Z());
 279 }
 280
 281 void FillPHOSHit(AliDetector* detector, AliHit* hit,
 282                  TH1F* hx, TH1F* hy, TH1F* hz)
 283 {
 284   // PHOS special
 285   Float_t id = ((AliPHOSHit*)hit)->GetId();
 286
 287   TVector3 pos;
 288   ((AliPHOS*)detector)->GetGeometry()->RelPosInAlice(id, pos);
 289
 290   if (hx) hx->Fill(pos.X());
 291   if (hy) hy->Fill(pos.Y());
 292   if (hz) hz->Fill(pos.Z());
 293 }
 294
 295 void DrawHistograms(TString& title, TString& gifName,
 296                     TH1F* h1, TH1F* h2, TH1F* h3, TH1F* h4)
 297 {
 298 // Create canvas, set the view range, show histograms.
 299 // ---
 300
 301   if (h1 && h2 && h3 && h4) {
 302     // create canvas
 303     TCanvas* canvas1 = new TCanvas("c1", title, 400, 10, 400, 300);
 304     TCanvas* canvas2 = new TCanvas("c2", title, 400, 10, 400, 300);
 305     TCanvas* canvas3 = new TCanvas("c3", title, 400, 10, 400, 300);
 306     TCanvas* canvas4 = new TCanvas("c4", title, 400, 10, 400, 300);
 307     //canvas->Divide(2,2);
 308
 309     // draw histograms
 310     canvas1->cd(1); h1->Draw();
 311     canvas2->cd(1); h2->Draw();
 312     canvas3->cd(1); h3->Draw();
 313     canvas4->cd(1); h4->Draw();
 314
 315     // save gif
 316     //canvas->SaveAs(gifNameBase + "_x.gif");
 317     canvas1->SaveAs(gifName+"1");
 318     canvas2->SaveAs(gifName+"2");
 319     canvas3->SaveAs(gifName+"3");
 320     canvas4->SaveAs(gifName+"4");
 321   }
 322 }
 323
 324 void DrawHistogram(TString& title, TString& gifName, TH1F* h)
 325 {
 326 // Create canvas, set the view range, show histograms.
 327 // ---
 328
 329   if (h){
 330     // create canvas
 331     TCanvas* canvas = new TCanvas("c", title, 400, 10, 400, 300);
 332
 333     // draw histograms
 334     h->Draw();
 335
 336     // save gif
 337     canvas->SaveAs(gifName);
 338   }
 339 }