MUON/MUONocc.C

   1 void MUONocc (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
  12    if (gClassTable->GetID("AliRun") < 0) {
  13       gROOT->LoadMacro("loadlibs.C");
  14       loadlibs();
  15    }
  16 // Connect the Root Galice file containing Geometry, Kine and Hits
  17
  18     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  19
  20
  21     if (!file) {
  22         printf("\n Creating galice.root \n");
  23         file = new TFile("galice.root");
  24     } else {
  25         printf("\n galice.root found in file list");
  26     }
  27     file->ls();
  28
  29 // Get AliRun object from file or create it if not on file
  30     if (!gAlice) {
  31         gAlice = (AliRun*)(file->Get("gAlice"));
  32         if (gAlice) printf("AliRun object found on file\n");
  33         if (!gAlice) {
  34             printf("\n create new gAlice object");
  35             gAlice = new AliRun("gAlice","Alice test program");
  36         }
  37     }
  38
  39 //  Create some histograms
  40
  41     TH1F *Hits[10], *HitDensity[10], *Occ0[10],*Occ1[10], *Mult[10];
  42     for (Int_t i=0; i<10; i++) {
  43         Hits[i]     =  new TH1F("hits1","Hit Distribution",26,0,260);
  44         Hits[i]->SetFillColor(0);
  45         Hits[i]->SetXTitle("R (cm)");
  46
  47         HitDensity[i]  =  new TH1F("dhits1","Hit Density Distribution",26,0,260);
  48         HitDensity[i]->SetFillColor(0);
  49         HitDensity[i]->SetXTitle("R (cm)");
  50
  51         Occ0[i]  =  new TH1F("occ0","Occupancy Density",26,0,260);
  52         Occ0[i] -> SetFillColor(0);
  53         Occ0[i] -> SetXTitle("R (cm)");
  54         Occ1[i]  =  new TH1F("occ1","Occupancy",26,0,260);
  55         Occ1[i] -> SetFillColor(0);
  56         Occ1[i] -> SetXTitle("R (cm)");
  57         Occ1[i] -> SetYTitle("Occupancy");
  58         Mult[i]  =  new TH1F("mult","Mult distribution",26,0,260);
  59         Mult[i] -> SetFillColor(0);
  60         Mult[i] -> SetXTitle("R (cm)");
  61     }
  62
  63
  64    TH1F *theta   =  new TH1F("theta","Theta distribution",180,0,180);
  65
  66    TH1F *emult   =  new TH1F("emult","Event Multiplicity",100,0,1000);
  67    AliMUONChamber*  iChamber;
  68    AliMUONSegmentation* seg;
  69
  70 //
  71 //   Loop over events
  72 //
  73    Float_t dpx[2], dpy[2];
  74    Int_t mdig =0;
  75    Int_t nhit=0;
  76
  77    for (Int_t nev=0; nev<= evNumber2; nev++) {
  78        Int_t nparticles = gAlice->GetEvent(nev);
  79        if (nev < evNumber1) continue;
  80        if (nparticles <= 0) return;
  81
  82        AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
  83        printf("\n track %d %d \n ", nev, MUON);
  84
  85        TTree *TH = gAlice->TreeH();
  86        Int_t ntracks = TH->GetEntries();
  87 //
  88 //   Loop over events
  89 //
  90        Int_t EvMult=0;
  91
  92        for (Int_t track=0; track<ntracks;track++) {
  93            gAlice->ResetHits();
  94            Int_t nbytes += TH->GetEvent(track);
  95            if (MUON)  {
  96 //
  97 //  Loop over hits
  98 //
  99                Int_t NperCh=0;
 100
 101                for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
 102                    mHit;
 103                    mHit=(AliMUONHit*)MUON->NextHit())
 104                {
 105                    Int_t   nch   = mHit->Chamber();  // chamber number
 106                    Float_t x     = mHit->X();        // x-pos of hit
 107                    Float_t y     = mHit->Y();        // y-pos
 108                    Float_t Eloss = mHit->Eloss();
 109                    Float_t Theta = mHit->Theta();
 110                    theta->Fill(Theta,(float) 1);
 111                    if (nch >10) continue;
 112                    if (nch ==1) EvMult++;
 113 //
 114 //
 115                    iChamber = & MUON->Chamber(nch-1);
 116                    response= iChamber.GetResponseModel();
 117                    seg     = iChamber.GetSegmentationModel(1);
 118                    NperCh++;
 119                    Int_t i,j;
 120                    seg->GetPadIxy(x,y,i,j);
 121                    Int_t isec = seg->Sector(i,j);
 122                    if (isec ==1 && nch==ic) nhit++;
 123
 124                    Float_t a=seg->Dpx(isec)*seg->Dpy(isec);
 125                    Float_t r=TMath::Sqrt(x*x+y*y);
 126                    Float_t wgt=1/(2*10*TMath::Pi()*r)/(evNumber2+1);
 127                    wgt=wgt*(1.+24./(2.*TMath::Pi()*r));
 128                    Hits[nch-1]->Fill(r,(float) 1);
 129                    HitDensity[nch-1]->Fill(r,wgt);
 130                    Occ0[nch-1]->Fill(r,wgt*a);
 131                } // hit loop
 132            } // if MUON
 133        } // track loop
 134        emult->Fill(Float_t(EvMult), (Float_t) 1);
 135
 136        Int_t iseg=1;
 137
 138        for (Int_t ich=0; ich<10; ich++) {
 139            iChamber = & MUON->Chamber(ich);
 140            seg=iChamber.GetSegmentationModel(iseg);
 141            gAlice->ResetDigits();
 142            gAlice->TreeD()->GetEvent(iseg);
 143            TClonesArray *MUONDigits  = MUON->DigitsAddress(ich);
 144            Int_t Ndigits=MUONDigits->GetEntriesFast();
 145            AliMUONDigit* dig;
 146            printf("\n Reading %d digits\n", Ndigits);
 147            if (MUONDigits)  {
 148                for (Int_t ndig=0; ndig<Ndigits; ndig++)
 149                {
 150                    dig = (AliMUONDigit*)MUONDigits->UncheckedAt(ndig);
 151                    Int_t i=dig->PadX();
 152                    Int_t j=dig->PadY();
 153                    Float_t x,y;
 154                    seg->GetPadCxy(i,j,x,y);
 155                    Int_t isec = seg->Sector(i,j);
 156                    Float_t a=seg->Dpx(isec)*seg->Dpy(isec);
 157                    Float_t r=TMath::Sqrt(x*x+y*y);
 158 //             if (r<25)
 159 //             printf("\n Sector,a %d %f", isec,a);
 160                    if (isec==1) mdig++;
 161
 162                    Float_t wgt;
 163                    wgt=1/(2.*10.*TMath::Pi()*r)/(evNumber2+1)*a;
 164 // Take into account inefficiency due to frames
 165                    wgt=wgt*(1.+24./(2.*TMath::Pi()*r));
 166                    Occ1[ich]->Fill(r,wgt);
 167                } // digit loop
 168                Mult[ich]->Divide(Occ1[ich],Occ0[ich]);
 169            } // chamber loop
 170        } // if MUONDigits
 171    } // event loop
 172 //
 173    printf("\n hits, digits %d %d\n ", nhit, mdig);
 174
 175
 176
 177 //Create a canvas, set the view range, show histograms
 178    TCanvas *c1 = new TCanvas("c1","Hit Densities",400,10,600,700);
 179    pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
 180    pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
 181    pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
 182    pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
 183    pad11->SetFillColor(0);
 184    pad12->SetFillColor(0);
 185    pad13->SetFillColor(0);
 186    pad14->SetFillColor(0);
 187    pad11->Draw();
 188    pad12->Draw();
 189    pad13->Draw();
 190    pad14->Draw();
 191
 192    pad11->cd();
 193    HitDensity[0]->Draw();
 194    pad12->cd();
 195    HitDensity[1]->Draw();
 196    pad13->cd();
 197    HitDensity[2]->Draw();
 198    pad14->cd();
 199    HitDensity[3]->Draw();
 200
 201    TCanvas *c2 = new TCanvas("c2","Hit Densities",400,10,600,700);
 202    pad21 = new TPad("pad21"," ",0.01,0.51,0.49,0.99);
 203    pad22 = new TPad("pad22"," ",0.51,0.51,0.99,0.99);
 204    pad23 = new TPad("pad23"," ",0.01,0.01,0.49,0.49);
 205    pad24 = new TPad("pad24"," ",0.51,0.01,0.99,0.49);
 206    pad21->SetFillColor(0);
 207    pad22->SetFillColor(0);
 208    pad23->SetFillColor(0);
 209    pad24->SetFillColor(0);
 210    pad21->Draw();
 211    pad22->Draw();
 212    pad23->Draw();
 213    pad24->Draw();
 214
 215    pad21->cd();
 216    HitDensity[4]->Draw();
 217    pad22->cd();
 218    HitDensity[5]->Draw();
 219    pad23->cd();
 220    HitDensity[6]->Draw();
 221    pad24->cd();
 222    HitDensity[7]->Draw();
 223
 224    TCanvas *c3 = new TCanvas("c3","Hit Densities",400,10,600,700);
 225    pad31 = new TPad("pad31"," ",0.01,0.51,0.49,0.99);
 226    pad32 = new TPad("pad32"," ",0.51,0.51,0.99,0.99);
 227    pad33 = new TPad("pad33"," ",0.01,0.01,0.49,0.49);
 228    pad34 = new TPad("pad34"," ",0.51,0.01,0.99,0.49);
 229    pad31->SetFillColor(0);
 230    pad32->SetFillColor(0);
 231    pad33->SetFillColor(0);
 232    pad34->SetFillColor(0);
 233    pad31->Draw();
 234    pad32->Draw();
 235    pad33->Draw();
 236    pad34->Draw();
 237
 238    pad31->cd();
 239    HitDensity[8]->Draw();
 240    pad32->cd();
 241    HitDensity[9]->Draw();
 242
 243    TCanvas *c4 = new TCanvas("c4","Occupancies",400,10,600,700);
 244    pad41 = new TPad("pad41"," ",0.01,0.51,0.49,0.99);
 245    pad42 = new TPad("pad42"," ",0.51,0.51,0.99,0.99);
 246    pad43 = new TPad("pad43"," ",0.01,0.01,0.49,0.49);
 247    pad44 = new TPad("pad44"," ",0.51,0.01,0.99,0.49);
 248    pad41->SetFillColor(0);
 249    pad42->SetFillColor(0);
 250    pad43->SetFillColor(0);
 251    pad44->SetFillColor(0);
 252    pad41->Draw();
 253    pad42->Draw();
 254    pad43->Draw();
 255    pad44->Draw();
 256
 257
 258    pad41->cd();
 259    Occ1[0]->Draw();
 260    pad42->cd();
 261    Occ1[2]->Draw();
 262    pad43->cd();
 263    Occ1[4]->Draw();
 264    pad44->cd();
 265    Occ1[6]->Scale(1.25);
 266    Occ1[6]->Draw();
 267
 268    TCanvas *c5 = new TCanvas("c5","Occupancies",400,10,600,700);
 269    pad51 = new TPad("pad41"," ",0.01,0.51,0.49,0.99);
 270    pad52 = new TPad("pad42"," ",0.51,0.51,0.99,0.99);
 271    pad53 = new TPad("pad43"," ",0.01,0.01,0.49,0.49);
 272    pad54 = new TPad("pad44"," ",0.51,0.01,0.99,0.49);
 273    pad51->SetFillColor(0);
 274    pad52->SetFillColor(0);
 275    pad53->SetFillColor(0);
 276    pad54->SetFillColor(0);
 277    pad51->Draw();
 278    pad52->Draw();
 279    pad53->Draw();
 280    pad54->Draw();
 281
 282
 283    pad51->cd();
 284    Occ1[8]->Scale(1.25);
 285
 286    Occ1[8]->Draw();
 287 }
 288
 289
 290
 291
 292
 293
 294
 295
 296
 297