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