AliGeant4/macro/hits.C

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