MUON/MUONhits.C

   1 void MUONhits (Int_t evNumber1=0,Int_t evNumber2=0, Int_t ic=1)
   2 {
   3 /////////////////////////////////////////////////////////////////////////
   4 //   This macro is a small example of a ROOT macro
   5 //   illustrating how to read the output of GALICE
   6 //   and do some analysis.
   7 //
   8 /////////////////////////////////////////////////////////////////////////
   9
  10 // Dynamically link some shared libs
  11     Float_t rmin[14] = {17.5, 17.5, 23.5, 23.5, 33.5, 33.5, 43., 43., 50., 50,
  12                         56.1, 56.1, 59.6, 59.6};
  13     Float_t rmax[14] = {81.6, 81.6, 109.3, 109.3, 154.4, 154.4, 197.8,
  14                         197.8, 229.5, 229.5, 254., 254, 270., 270.};
  15
  16
  17     if (gClassTable->GetID("AliRun") < 0) {
  18         gROOT->LoadMacro("loadlibs.C");
  19         loadlibs();
  20     }
  21 // Connect the Root Galice file containing Geometry, Kine and Hits
  22
  23     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  24
  25
  26     if (!file) {
  27         printf("\n Creating galice.root \n");
  28         file = new TFile("galice.root");
  29     } else {
  30         printf("\n galice.root found in file list");
  31     }
  32     file->ls();
  33 //
  34     TDatabasePDG*  DataBase = new TDatabasePDG();
  35
  36 // Get AliRun object from file or create it if not on file
  37     if (!gAlice) {
  38         gAlice = (AliRun*)(file->Get("gAlice"));
  39         if (gAlice) printf("AliRun object found on file\n");
  40         if (!gAlice) {
  41             printf("\n create new gAlice object");
  42             gAlice = new AliRun("gAlice","Alice test program");
  43         }
  44     }
  45
  46 //  Create some histograms
  47
  48     TH1F *Hits[14],   *HitDensity[14];
  49     TH1F *Hits_g[14], *HitDensity_g[14];
  50     TH1F *Hits_n[14], *HitDensity_n[14];
  51     TH1F *Mom[14], *Mom_g[14], *Mom_n[14];
  52
  53     for (Int_t i=0; i<14; i++) {
  54         Hits[i]     =  new TH1F("hits1","Hit Distribution",30,0,300);
  55         Hits[i]->SetFillColor(0);
  56         Hits[i]->SetXTitle("R (cm)");
  57
  58         HitDensity[i]  =  new TH1F("dhits1","Hit Density",30,0,300);
  59         HitDensity[i]->SetFillColor(0);
  60         HitDensity[i]->SetXTitle("R (cm)");
  61
  62         Hits_g[i]     =  new TH1F("hits1","Hit Distribution",30,0,300);
  63         Hits_g[i]->SetFillColor(0);
  64         Hits_g[i]->SetXTitle("R (cm)");
  65
  66         HitDensity_g[i]  =  new TH1F("dhits1","Hit Density",30,0,300);
  67         HitDensity_g[i]->SetFillColor(0);
  68         HitDensity_g[i]->SetXTitle("R (cm)");
  69
  70         Mom[i]  =  new TH1F("mom","Energy",70,-6,1);
  71         Mom[i]->SetFillColor(0);
  72         Mom[i]->SetXTitle("E (GeV)");
  73
  74         Mom_g[i]  =  new TH1F("mom","Energy",70,-6,1);
  75         Mom_g[i]->SetFillColor(0);
  76         Mom_g[i]->SetXTitle("E (GeV)");
  77
  78         Mom_n[i]  =  new TH1F("mom","Energy",70,-6,1);
  79         Mom_n[i]->SetFillColor(0);
  80         Mom_n[i]->SetXTitle("E (GeV)");
  81
  82         Hits_n[i]     =  new TH1F("hits1","Hit Distribution",30,0,300);
  83         Hits_n[i]->SetFillColor(0);
  84         Hits_n[i]->SetXTitle("R (cm)");
  85
  86         HitDensity_n[i]  =  new TH1F("dhits1","Hit Density",30,0,300);
  87         HitDensity_n[i]->SetFillColor(0);
  88         HitDensity_n[i]->SetXTitle("R (cm)");
  89     }
  90
  91     TH1F *theta   =  new TH1F("theta","Theta distribution",180,0,180);
  92
  93     TH1F *emult   =  new TH1F("emult","Event Multiplicity",100,0,1000);
  94     AliMUONChamber*  iChamber;
  95     AliSegmentation* seg;
  96
  97 //
  98 //   Loop over events
  99 //
 100     Float_t dpx[2], dpy[2];
 101     Int_t mdig =0;
 102     Int_t nhit=0;
 103
 104     for (Int_t nev=0; nev<= evNumber2; nev++) {
 105         Int_t nparticles = gAlice->GetEvent(nev);
 106         if (nev < evNumber1) continue;
 107         if (nparticles <= 0) return;
 108
 109         AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
 110         printf("\n track %d %d \n ", nev, MUON);
 111
 112         TTree *TH = gAlice->TreeH();
 113         Int_t ntracks = TH->GetEntries();
 114 //
 115 //   Loop over events
 116 //
 117         Int_t EvMult=0;
 118
 119         for (Int_t track=0; track<ntracks;track++) {
 120             gAlice->ResetHits();
 121             Int_t nbytes += TH->GetEvent(track);
 122             if (MUON)  {
 123 //
 124 //  Loop over hits
 125 //
 126                 Int_t NperCh=0;
 127
 128                 for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
 129                     mHit;
 130                     mHit=(AliMUONHit*)MUON->NextHit())
 131                 {
 132                     Int_t   nch   =    mHit->Chamber();  // chamber number
 133                     Float_t x     =    mHit->X();        // x-pos of hit
 134                     Float_t y     =    mHit->Y();        // y-pos
 135                     Float_t Eloss =    mHit->Eloss();
 136                     Float_t Theta =    mHit->Theta();
 137                     Float_t Particle = mHit->Particle();
 138                     Float_t P     =    mHit->Momentum();
 139                     TParticlePDG* Part = DataBase->GetParticle(Particle);
 140                     Double_t mass = Part->Mass();
 141
 142                     if (nch >13) continue;
 143                     if (nch ==1) EvMult++;
 144                     if (mHit->Age() > 5.e-6) continue;
 145
 146                     Float_t r=TMath::Sqrt(x*x+y*y);
 147                     if (nch ==0) continue;
 148
 149                     if (r < rmin[nch-1]) continue;
 150                     if (r > rmax[nch-1]) continue;
 151
 152                     Float_t wgt=1/(2*10*TMath::Pi()*r)/(evNumber2+1);
 153
 154                     Float_t Ekin=TMath::Sqrt(P*P+mass*mass)-mass;
 155                     if (Particle == 2112) {
 156                         Hits_n[nch-1]->Fill(r,(float) 1);
 157                         HitDensity_n[nch-1]->Fill(r,wgt);
 158                         Mom_n[nch-1]->Fill(TMath::Log10(Ekin),1);
 159                     } else if (Particle == 22) {
 160                         Hits_g[nch-1]->Fill(r,(float) 1);
 161                         HitDensity_g[nch-1]->Fill(r,wgt);
 162                         Mom_g[nch-1]->Fill(TMath::Log10(Ekin),1);
 163                     } else {
 164                         Hits[nch-1]->Fill(r,(float) 1);
 165                         HitDensity[nch-1]->Fill(r,wgt);
 166                         Mom[nch-1]->Fill(TMath::Log10(Ekin),1);
 167                     }
 168                 } // hit loop
 169             } // if MUON
 170         } // track loop
 171     }
 172
 173 //Create a canvas, set the view range, show histograms
 174     Int_t k;
 175     TCanvas *c1 = new TCanvas("c1","Hit Densities",400,10,600,700);
 176     c1->Divide(2,4);
 177     for (k=0; k<7; k++) {
 178         c1->cd(k+1);
 179         HitDensity[2*k]->Draw();
 180     }
 181
 182     TCanvas *c2 = new TCanvas("c2","Hit Densities (gamma)",400,10,600,700);
 183     c2->Divide(2,4);
 184     for (k=0; k<7; k++) {
 185         c2->cd(k+1);
 186         HitDensity_g[2*k]->Draw();
 187     }
 188
 189     TCanvas *c3 = new TCanvas("c3","Hit Densities (neutron)",400,10,600,700);
 190     c3->Divide(2,4);
 191     for (k=0; k<7; k++) {
 192         c3->cd(k+1);
 193         HitDensity_n[2*k]->Draw();
 194     }
 195
 196     TCanvas *c4 = new TCanvas("c4","Energy (charged)",400,10,600,700);
 197     c4->Divide(2,4);
 198     for (k=0; k<7; k++) {
 199         c4->cd(k+1);
 200         Mom[2*k]->Draw();
 201     }
 202
 203     TCanvas *c5 = new TCanvas("c5","Energy (gamma)",400,10,600,700);
 204     c5->Divide(2,4);
 205     for (k=0; k<7; k++) {
 206         c5->cd(k+1);
 207         Mom_g[2*k]->Draw();
 208     }
 209
 210     TCanvas *c6 = new TCanvas("c6","Energy (gamma)",400,10,600,700);
 211     c6->Divide(2,4);
 212     for (k=0; k<7; k++) {
 213         c6->cd(k+1);
 214         Mom_n[2*k]->Draw();
 215     }
 216 }
 217
 218
 219
 220
 221
 222
 223