MUON/MUONacc.C

   1 Bool_t SelectThis(AliMUONHit* mHit);
   2
   3 Bool_t SelectThis(Float_t x, Float_t y, Float_t d)
   4 {
   5 //    Float_t x     =    mHit->fX;        // x-pos
   6 //    Float_t y     =    mHit->fY;        // y-pos
   7     if (TMath::Abs(x) > d/2. && TMath::Abs(y) > d/2.) {
   8         return 1;
   9     } else {
  10         return 0;
  11     }
  12 }
  13
  14
  15 void MUONacc (Float_t d=3., Int_t evNumber1=0, Int_t evNumber2=0)
  16 {
  17 // Dynamically link some shared libs
  18     Float_t rmin[14] = {17.5, 17.5, 23.5, 23.5, 33.5, 33.5, 43., 43., 50., 50,
  19                         56.1, 56.1, 59.6, 59.6};
  20     Float_t rmax[14] = {91.5, 91.5, 122.5, 122.5, 158.3, 158.3, 260.0,
  21                         260.0, 260.0, 260.0, 293., 293, 293., 293.};
  22
  23
  24     if (gClassTable->GetID("AliRun") < 0) {
  25         gROOT->LoadMacro("loadlibs.C");
  26         loadlibs();
  27     }
  28 // Connect the Root Galice file containing Geometry, Kine and Hits
  29
  30     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  31
  32
  33     if (!file) {
  34         printf("\n Creating galice.root \n");
  35         file = new TFile("galice.root");
  36     } else {
  37         printf("\n galice.root found in file list");
  38     }
  39     file->ls();
  40 //
  41     TDatabasePDG*  DataBase = new TDatabasePDG();
  42
  43 // Get AliRun object from file or create it if not on file
  44     if (!gAlice) {
  45         gAlice = (AliRun*)(file->Get("gAlice"));
  46         if (gAlice) printf("AliRun object found on file\n");
  47         if (!gAlice) {
  48             printf("\n create new gAlice object");
  49             gAlice = new AliRun("gAlice","Alice test program");
  50         }
  51     }
  52
  53 //  Create some histograms
  54
  55     TH1F *theta   =  new TH1F("theta","Theta distribution",180,0,180);
  56     TH1F *emult   =  new TH1F("emult","Event Multiplicity",100,0,1000);
  57
  58     AliMUONChamber*  iChamber;
  59     AliSegmentation* seg;
  60
  61 //
  62 //   Loop over events
  63 //
  64
  65     for (Int_t nev=0; nev<= evNumber2; nev++) {
  66         Int_t nparticles = gAlice->GetEvent(nev);
  67         if (nev < evNumber1) continue;
  68         if (nparticles <= 0) return;
  69
  70         AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
  71
  72
  73         TTree *TH = gAlice->TreeH();
  74         Int_t ntracks = TH->GetEntries();
  75         TClonesArray *fPartArray = gAlice->Particles();
  76         Int_t npart = fPartArray->GetEntriesFast();
  77
  78         Int_t EvMult=0;
  79         Int_t Nsel=0;
  80 //
  81 // Loop over primary particles (jpsi. upsilon, ...)
  82 //
  83         Int_t prim=0;
  84
  85         for (Int_t part=0; part<3000; part+=3) {
  86             TParticle *MPart = (TParticle*) fPartArray->UncheckedAt(part);
  87             Int_t mpart = MPart->GetPdgCode();
  88             Int_t child1 = MPart->GetFirstDaughter();
  89             Int_t child2 = MPart->GetLastDaughter();
  90             Int_t mother = MPart->GetFirstMother();
  91 //
  92 // Loop over muons
  93 //
  94             Int_t muon1=2*prim;
  95             Int_t muon2=muon1+1;
  96
  97 //          printf("\n PDG Code %d %d %d %d %d %d %d %d \n",
  98 //                 part, mpart, child1, child2, mother,prim, muon1,muon2);
  99
 100
 101             Bool_t selected[2]={1,1};
 102             Int_t muons=0;
 103             for (Int_t track=muon1; track<=muon2;track++) {
 104                 gAlice->ResetHits();
 105                 Int_t nbytes += TH->GetEvent(track);
 106
 107                 if (MUON)  {
 108 //
 109 //  Loop over hits
 110 //
 111                     for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
 112                         mHit;
 113                         mHit=(AliMUONHit*)MUON->NextHit())
 114                     {
 115                         Int_t   nch   =    mHit->fChamber;  // chamber number
 116                         Float_t x     =    mHit->fX;        // x-pos of hit
 117                         Float_t y     =    mHit->fY;        // y-pos
 118                         Float_t z     =    mHit->fZ;        // z-pos
 119                         Float_t Eloss =    mHit->fEloss;    // energy loss
 120                         Float_t Theta =    mHit->fTheta;    // theta
 121                         Float_t Particle = mHit->fParticle; // Particle type
 122                         Int_t itrack   = Int_t(mHit->fTrack);
 123                         TParticle *thePart = (TParticle*) fPartArray->UncheckedAt(itrack);
 124                         Float_t pTheta=thePart->Theta();
 125
 126                         if (Particle != kMuonPlus && Particle != kMuonMinus) continue;
 127
 128                         Float_t P    =
 129                             TMath::Sqrt(mHit->fCxHit*mHit->fCxHit+
 130                                         mHit->fCyHit*mHit->fCyHit+
 131                                         mHit->fCzHit*mHit->fCzHit);
 132                         Float_t R   = TMath::Sqrt(x*x+y*y);
 133                         TParticlePDG* Part = DataBase->GetParticle(Particle);
 134                         Double_t mass = Part->Mass();
 135                         theta->Fill(pTheta*180./TMath::Pi(),1.);
 136
 137                         if (nch>14) continue;
 138 //                      printf("\n %f %f %f %d\n ", R, rmin[nch-1], rmax[nch-1], nch);
 139 // Geometrical acceptance of frames
 140 //                      selected[muons] = selected[muons]&&SelectThis(x,y,d);
 141 // Cut on default rmin and rmax
 142 //                      selected[muons] = selected[muons] &&( R > rmin[nch-1] && R < rmax[nch-1]);
 143 // Cut on rmin
 144 //                      selected[muons] = selected[muons] &&( R > rmin[nch-1]);
 145 // Cut on rmax
 146 //                      selected[muons] = selected[muons] &&( R < rmax[nch-1]);
 147 //                      selected[muons] = selected[muons] &&( R > z*TMath::Tan(2.0*TMath::Pi()/180.));
 148 //                      Bool_t ok = (z<970 && ( R > z*TMath::Tan(2.0*TMath::Pi()/180.))) ||
 149 //                          (z>970 && ( R > 970 * TMath::Tan(2.0*TMath::Pi()/180.)
 150 //                                      +(z-970)* TMath::Tan(1.6*TMath::Pi()/180.)));
 151 //
 152                         Bool_t ok = (z<970 && ( R < z*TMath::Tan(9.0*TMath::Pi()/180.))) ||
 153                             (z>970 && ( R < 970 * TMath::Tan(9.0*TMath::Pi()/180.)
 154                                         +(z-970)* TMath::Tan(12.0*TMath::Pi()/180.)));
 155
 156
 157                         selected[muons] = selected[muons] && ok;
 158
 159                     } // hit loop
 160                 } // if MUON
 161                 muons++;
 162             } // track loop
 163             if (selected[0] && selected[1]) Nsel++;
 164             prim++;
 165         } // primary loop
 166         printf("\n Selected %d %d \n", Nsel, prim);
 167     }
 168
 169 //Create a canvas, set the view range, show histograms
 170     Int_t k;
 171     TCanvas *c1 = new TCanvas("c1","Canvas 1",400,10,600,700);
 172     c1->Divide(2,4);
 173     c1->cd(1);
 174     theta->Draw();
 175
 176 }
 177
 178
 179
 180
 181
 182
 183
 184
 185