MUON/MUONdoubles.C

   1 #include "iostream.h"
   2
   3 void MUONdoubles (Int_t evNumber1=0,Int_t evNumber2=0)
   4 {
   5     {
   6     hprof  = new TProfile("hprof","Profile dmin vs r",24,0,120,0,50);
   7     }
   8     hdist1  = new TH1F("hdist1","distance",100,0,10);
   9     hdist2  = new TH1F("hdist2","distance",100,0,10);
  10     Float_t a1=3.7/TMath::Sqrt(4);
  11     Float_t a2=26/TMath::Sqrt(4);
  12     for (Int_t i=1; i<10000; i++) {
  13         Float_t x1,x2,y1,y2;
  14         x1 = (2*gRandom->Rndm(i)-1.)*a1;
  15         x2 = (2*gRandom->Rndm(i)-1.)*a1;
  16         y1 = (2*gRandom->Rndm(i)-1.)*a1;
  17         y2 = (2*gRandom->Rndm(i)-1.)*a1;
  18         Float_t d=TMath::Sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
  19         hdist1->Fill(d,1);
  20     }
  21     Float_t xh[100], yh[100];
  22     for (Int_t k=0; k<1000; k++) {
  23         for (Int_t i=0; i<100; i++) {
  24             xh[i] = (2*gRandom->Rndm(i)-1.)*a2;
  25             yh[i] = (2*gRandom->Rndm(i)-1.)*a2;
  26         }
  27         Float_t x1,y1;
  28
  29         for (Int_t i=0; i<10; i++) {
  30             Float_t dmin=1000;
  31             x1 = (2*gRandom->Rndm(i)-1.)*a2;
  32             y1 = (2*gRandom->Rndm(i)-1.)*a2;
  33             for (Int_t j=0; j<100; j++) {
  34                 Float_t d=TMath::Sqrt((x1-xh[j])*(x1-xh[j])+(y1-yh[j])*(y1-yh[j]));
  35                 if (d<dmin) dmin=d;
  36             }
  37             hdist2->Fill(dmin,1);
  38         }
  39     }
  40
  41
  42
  43 /////////////////////////////////////////////////////////////////////////
  44 //   This macro is a small example of a ROOT macro
  45 //   illustrating how to read the output of GALICE
  46 //   and do some analysis.
  47 //
  48 /////////////////////////////////////////////////////////////////////////
  49
  50 // Dynamically link some shared libs
  51
  52     if (gClassTable->GetID("AliRun") < 0) {
  53         gROOT->LoadMacro("loadlibs.C");
  54         loadlibs();
  55     }
  56
  57 // Connect the Root Galice file containing Geometry, Kine and Hits
  58
  59     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  60     if (!file) file = new TFile("galice.root","UPDATE");
  61
  62 // Get AliRun object from file or create it if not on file
  63
  64     if (!gAlice) {
  65         gAlice = (AliRun*)file->Get("gAlice");
  66         if (gAlice) printf("AliRun object found on file\n");
  67         if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  68     }
  69 //
  70 // Set reconstruction models
  71 //
  72 // Get pointers to Alice detectors and Digits containers
  73     AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
  74 //
  75 // Book some histograms
  76     TH1F *Hcentre[10], *Hall[10], *Hcfrac[10], *Htail[10], *Htfrac[10];
  77     TH1F *Hdcut[10], *Hdall[10], *Hdfrac[10];
  78
  79     for (Int_t i=0; i<10; i++) {
  80         Hcentre[i]     =  new TH1F("Hcentre","Hit Distribution",26,0,260);
  81         Hcentre[i]->SetFillColor(0);
  82         Hcentre[i]->SetXTitle("R (cm)");
  83
  84         Hall[i]     =  new TH1F("Hall","Hit Distribution",26,0,260);
  85         Hall[i]->SetFillColor(0);
  86         Hall[i]->SetXTitle("R (cm)");
  87
  88         Hcfrac[i]     =  new TH1F("Hcfrag","Hit Distribution",26,0,260);
  89         Hcfrac[i]->SetFillColor(0);
  90         Hcfrac[i]->SetXTitle("R (cm)");
  91
  92         Htail[i]     =  new TH1F("Htail","Hit Distribution",26,0,260);
  93         Htail[i]->SetFillColor(0);
  94         Htail[i]->SetXTitle("R (cm)");
  95
  96         Htfrac[i]     =  new TH1F("Htfrag","Hit Distribution",26,0,260);
  97         Htfrac[i]->SetFillColor(0);
  98         Htfrac[i]->SetXTitle("R (cm)");
  99
 100         Hdall[i]     =  new TH1F("Hdall","Hit Distribution",26,0,260);
 101         Hdall[i]->SetFillColor(0);
 102         Hdall[i]->SetXTitle("R (cm)");
 103
 104         Hdcut[i]     =  new TH1F("Hdcut","Hit Distribution",26,0,260);
 105         Hdcut[i]->SetFillColor(0);
 106         Hdcut[i]->SetXTitle("R (cm)");
 107
 108         Hdfrac[i]     =  new TH1F("Hdfrac","Hit Distribution",26,0,260);
 109         Hdfrac[i]->SetFillColor(0);
 110         Hdfrac[i]->SetXTitle("R (cm)");
 111
 112     }
 113
 114 //
 115 //   Loop over events
 116 //
 117     Int_t Nh=0;
 118     Int_t Nh1=0;
 119     for (int nev=0; nev<= evNumber2; nev++) {
 120         Int_t nparticles = gAlice->GetEvent(nev);
 121         cout << "nev         " << nev <<endl;
 122         cout << "nparticles  " << nparticles <<endl;
 123         if (nev < evNumber1) continue;
 124         if (nparticles <= 0) return;
 125
 126         TTree *TH = gAlice->TreeH();
 127         Int_t ntracks = TH->GetEntries();
 128         cout<<"ntracks "<<ntracks<<endl;
 129
 130         Int_t nbytes = 0;
 131
 132
 133         TClonesArray *Particles = gAlice->Particles();
 134         TTree *TD = gAlice->TreeD();
 135         Int_t nent=TD->GetEntries();
 136         printf("Found %d entries in the tree (must be one per cathode per event!)\n",nent);
 137         Float_t x[2000];
 138         Float_t y[2000];
 139         Int_t nhit=0;
 140
 141         if (MUON) {
 142 //
 143 // Loop on chambers and on cathode planes
 144 //
 145             TTree *TH = gAlice->TreeH();
 146             Int_t ntracks = TH->GetEntries();
 147
 148 //
 149 //              Loop over Hits
 150 //
 151             for (i=0; i<1000; i++) {
 152                 Float_t r  =100.*gRandom->Rndm(i);
 153                 Float_t phi=2.*TMath::Pi()*gRandom->Rndm(i);
 154                 Float_t xr=r*TMath::Sin(phi);
 155                 Float_t yr=r*TMath::Cos(phi);
 156                 Hdall[0]->Fill(r,1.);
 157                 Float_t dmin=100000.;
 158
 159                 if (i==0) {
 160                     for (Int_t track=0; track<ntracks; track++) {
 161                         gAlice->ResetHits();
 162                         Int_t nbytes += TH->GetEvent(track);
 163                         for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
 164                             mHit;
 165                             mHit=(AliMUONHit*)MUON->NextHit())
 166                         {
 167                             Int_t   nch   = mHit->fChamber;  // chamber number
 168                             if (nch!=1) continue;
 169
 170                             x[nhit]     = mHit->fX;        // x-pos of hit
 171                             y[nhit]     = mHit->fY;        // y-pos
 172                             nhit++;
 173                         } // hit loop
 174                     } // track loop
 175                 } else {
 176                     for (Int_t nh=0; nh<nhit; nh++) {
 177                         Float_t d=TMath::Sqrt((x[nh]-xr)*(x[nh]-xr)+(y[nh]-yr)*(y[nh]-yr));
 178 //                      printf ("\n r,d %f %f",r,d);
 179                         Float_t dx=TMath::Abs(x[nh]-xr);
 180                         Float_t dy=TMath::Abs(y[nh]-yr);
 181                         if (d < dmin) dmin=d;
 182
 183                         if (dx<0.5 && dy < 0.75) Hdcut[0]->Fill(r,1.);
 184                     } // hit loop
 185                 } // i==0
 186                 if (r>20) hprof->Fill(r,dmin,1);
 187             } // random loop
 188             Int_t icat=0;
 189             gAlice->ResetDigits();
 190             gAlice->TreeD()->GetEvent(icat+1); // spurious +1 ...
 191
 192             for (Int_t ich=0;ich<1;ich++) {
 193                 AliMUONChamber iChamber= MUON->Chamber(ich);
 194                 TClonesArray *MUONdigits  = MUON->DigitsAddress(ich);
 195                 if (MUONdigits == 0) continue;
 196                 //
 197                 // Get ready the current chamber stuff
 198                 //
 199                 AliMUONResponse* response = iChamber.GetResponseModel();
 200                 AliMUONSegmentation*  seg = iChamber.GetSegmentationModel(icat+1);
 201                 HitMap = new AliMUONHitMapA1(seg, MUONdigits);
 202                 HitMap->FillHits();
 203                 Int_t nxmax=seg->Npx();
 204                 Int_t nymax=seg->Npy();
 205 //
 206 // generate random positions on the chamber
 207 //
 208                 for (Int_t i=0; i<4000; i++) {
 209                     Int_t ix = 2*Int_t((gRandom->Rndm(i)-0.5)*nxmax);
 210                     Int_t iy = 2*Int_t((gRandom->Rndm(i)-0.5)*nymax);
 211
 212 // test for centre overlap
 213 //
 214                     Float_t xp,yp;
 215                     seg->GetPadCxy(ix,iy,xp,yp);
 216                     Float_t r=TMath::Sqrt(xp*xp+yp*yp);
 217                     if (TMath::Abs(xp)>3 && TMath::Abs(yp)>3)
 218                         Hall[ich]->Fill(r,1.);
 219                     if (HitMap->TestHit(ix,iy) != 0) {
 220                         Hcentre[ich]->Fill(r,1.);
 221                     }
 222 // test for neighbour overlap
 223 //
 224                     Int_t nn;
 225                     Int_t X[20], Y[20];
 226                     seg->Neighbours(ix,iy,&nn,X,Y);
 227                     Bool_t hit=false;
 228
 229                     for (Int_t j=0; j<nn; j++) {
 230                         if (HitMap->TestHit(X[j],Y[j]) != 0) hit=true;
 231                     }
 232                     if (hit &&  HitMap->TestHit(ix,iy) == 0) Htail[ich]->Fill(r,1.);
 233                 } //random loop
 234                 delete HitMap;
 235             } // chamber loop
 236         }   // event loop
 237     } // if MUON
 238
 239     for (Int_t ich=0; ich<10; ich++) {
 240         Hcfrac[ich]->Divide(Hcentre[ich],Hall[ich]);
 241         Htfrac[ich]->Divide(Htail[ich],Hall[ich]);
 242         Hdfrac[ich]->Divide(Hdcut[ich],Hdall[ich]);
 243     }
 244
 245     TCanvas *c1 = new TCanvas("c1","Hit Densities",400,10,600,700);
 246     pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
 247     pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
 248     pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
 249     pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
 250     pad11->SetFillColor(0);
 251     pad12->SetFillColor(0);
 252     pad13->SetFillColor(0);
 253     pad14->SetFillColor(0);
 254     pad11->Draw();
 255     pad12->Draw();
 256     pad13->Draw();
 257     pad14->Draw();
 258
 259     pad11->cd();
 260     Hcentre[0]->Draw();
 261
 262     pad12->cd();
 263     Hall[0]->Draw();
 264
 265    pad13->cd();
 266    Hcfrac[0]->Draw();
 267
 268    pad14->cd();
 269    Htfrac[0]->Draw();
 270
 271    TCanvas *c2 = new TCanvas("c2","Hit Densities",400,10,600,700);
 272    pad21 = new TPad("pad21"," ",0.01,0.51,0.49,0.99);
 273    pad22 = new TPad("pad22"," ",0.51,0.51,0.99,0.99);
 274    pad23 = new TPad("pad23"," ",0.01,0.01,0.49,0.49);
 275    pad24 = new TPad("pad24"," ",0.51,0.01,0.99,0.49);
 276    pad21->SetFillColor(0);
 277    pad22->SetFillColor(0);
 278    pad23->SetFillColor(0);
 279    pad24->SetFillColor(0);
 280    pad21->Draw();
 281    pad22->Draw();
 282    pad23->Draw();
 283    pad24->Draw();
 284
 285    pad21->cd();
 286    hdist1->Draw();
 287
 288    pad22->cd();
 289    hdist2->Draw();
 290
 291    pad23->cd();
 292    Hdfrac[0]->Draw();
 293
 294    pad24->cd();
 295    hprof->Draw();
 296    file->Close();
 297 }
 298