MUON/MUONtest.C

   1 void MUONtest (Int_t evNumber1=0,Int_t evNumber2=0)
   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         gSystem->Load("$ALITOP/cern.so/lib/libpdfDUMMY.so");
  14         gSystem->Load("$ALITOP/cern.so/lib/libPythia.so");
  15         gSystem->Load("$ROOTSYS/lib/libEG.so");
  16         gSystem->Load("$ROOTSYS/lib/libEGPythia.so");
  17         gSystem->Load("libGeant3Dummy.so");    //a dummy version of Geant3
  18         gSystem->Load("PHOS/libPHOSdummy.so"); //the standard Alice classes
  19         gSystem->Load("libgalice.so");         // the standard Alice classes
  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 // Get AliRun object from file or create it if not on file
  35     if (!gAlice) {
  36         gAlice = (AliRun*)(file->Get("gAlice"));
  37         if (gAlice) printf("AliRun object found on file\n");
  38         if (!gAlice) {
  39             printf("\n create new gAlice object");
  40             gAlice = new AliRun("gAlice","Alice test program");
  41         }
  42     }
  43
  44 //  Create some histograms
  45
  46
  47    TH1F *ccharge1 = new TH1F("ccharge1","Cluster Charge distribution"
  48                            ,100,0.,500.);
  49    TH1F *pcharge1 = new TH1F("pcharge1","Pad      Charge distribution"
  50                            ,100,0.,200.);
  51    TH1F *xresid1  = new TH1F("xresid1","x-Residuals"
  52                            ,100,-2.0,2.0);
  53    TH1F *yresid1  = new TH1F("yresid1","y-Residuals"
  54                            ,100,-0.2,0.2);
  55    TH1F *npads1   = new TH1F("npads1" ,"Pads per Hit"
  56                            ,20,-0.5,19.5);
  57    TH2F *xresys1  = new TH2F("xresys1","x-Residuals systematics"
  58                            ,50,-0.4,0.4,100,-0.2,0.2);
  59    TH2F *yresys1  = new TH2F("yresys1","y-Residuals systematics"
  60                            ,50,-0.4,0.4,100,-0.1,0.1);
  61
  62    TH1F *ccharge2 = new TH1F("ccharge2","Cluster Charge distribution"
  63                            ,100,0.,500.);
  64    TH1F *pcharge2 = new TH1F("pcharge2","Pad      Charge distribution"
  65                            ,100,0.,200.);
  66    TH1F *xresid2  = new TH1F("xresid2","x-Residuals"
  67                            ,100,-2.0,2.0);
  68    TH1F *yresid2  = new TH1F("yresid2","y-Residuals"
  69                            ,100,-0.2,0.2);
  70    TH1F *npads2   = new TH1F("npads2" ,"Pads per Hit"
  71                            ,20,-0.5,19.5);
  72    TH2F *xresys2  = new TH2F("xresys2","x-Residuals systematics"
  73                            ,50,-0.4,0.4,100,-0.2,0.2);
  74    TH2F *yresys2  = new TH2F("yresys2","y-Residuals systematics"
  75                            ,50,-0.4,0.4,100,-0.1,0.1);
  76
  77
  78    AliMUONchamber*  iChamber;
  79 //
  80 //   Loop over events
  81 //
  82    Int_t Nh=0;
  83    Int_t Nh1=0;
  84    for (Int_t nev=0; nev<= evNumber2; nev++) {
  85        cout << "nev         " << nev <<endl;
  86        Int_t nparticles = gAlice->GetEvent(nev);
  87        cout << "nparticles  " << nparticles <<endl;
  88        if (nev < evNumber1) continue;
  89        if (nparticles <= 0) return;
  90
  91        AliMUON *MUON  = (AliMUON*) gAlice->GetDetector("MUON");
  92            printf("\n track %d %d \n ", nev, MUON);
  93
  94        TTree *TH = gAlice->TreeH();
  95        Int_t ntracks = TH->GetEntries();
  96        Int_t Nc=0;
  97        Int_t   npad[2];
  98        Float_t Q[2], xmean[2],ymean[2],xres[2],yres[2], xonpad[2], yonpad[2];
  99
 100 //
 101 //   Loop over events
 102 //
 103        for (Int_t track=0; track<ntracks;track++) {
 104
 105            gAlice->ResetHits();
 106            Int_t nbytes += TH->GetEvent(track);
 107            if (MUON)  {
 108                for(AliMUONhit* mHit=(AliMUONhit*)MUON->FirstHit(-1);
 109                    mHit;
 110                    mHit=(AliMUONhit*)MUON->NextHit())
 111                {
 112                    Int_t   nch   = mHit->fChamber;  // chamber number
 113                    Float_t x     = mHit->fX;        // x-pos of hit
 114                    Float_t y     = mHit->fY;        // y-pos
 115 //
 116 //
 117                    iChamber = & MUON->Chamber(nch-1);
 118                    response=iChamber->GetResponseModel();
 119 //
 120 //
 121                    if (nch <= 1) {
 122                        for (Int_t i=0; i<2; i++) {
 123                            xmean[i]=0;
 124                            ymean[i]=0;
 125                            Q[i]=0;
 126                            npad[i]=0;
 127                        }
 128                        for (AliMUONcluster *mPad=(AliMUONcluster*)MUON->FirstPad(mHit);
 129                             mPad;
 130                             mPad=(AliMUONcluster*)MUON->NextPad()
 131                            )
 132                        {
 133                            Int_t nseg     = mPad->fCathode;   // segmentation
 134                            Int_t nhit     = mPad->fHitNumber; // hit number
 135                            Int_t qtot     = mPad->fQ;         // charge
 136                            Int_t ipx      = mPad->fPadX;      // pad number on X
 137                            Int_t ipy      = mPad->fPadY;      // pad number on Y
 138                            Int_t iqpad    = mPad->fQpad;      // charge per pad
 139                            Int_t rpad     = mPad->fRpad;      // r-pos of pad
 140 //
 141 //
 142                            segmentation=iChamber->GetSegmentationModel(nseg);
 143                            Int_t ind=nseg-1;
 144 //                     printf("Pad of hit, padx, pady, iqpad, ncha %d %d %d %d %d %d\n",
 145 //                            nhit, ipx, ipy, iqpad, nseg, ind);
 146
 147
 148                            if (nseg == 1) {
 149                                pcharge1->Fill(iqpad,(float) 1);
 150                                ccharge1->Fill(qtot ,(float) 1);
 151                            }  else {
 152                                pcharge2->Fill(iqpad,(float) 1);
 153                                ccharge2->Fill(qtot ,(float) 1);
 154                            }
 155 // Calculate centre of gravity
 156 //
 157                            if (iqpad > 0) {
 158                                npad[ind]++;
 159                                xmean[ind]+=iqpad*segmentation->GetPadx(ipx);
 160                                ymean[ind]+=iqpad*segmentation->GetPady(ipy);
 161                                Q[ind]+=iqpad;
 162                            }
 163
 164                        } //pad hit loop
 165                        for (Int_t i=0; i<2; i++) {
 166                            segmentation = iChamber->GetSegmentationModel(i+1);
 167                            if (Q[i] >0) {
 168                                xmean[i] = xmean[i]/Q[i];
 169                                xres[i]  = xmean[i]-x;
 170                                ymean[i] = ymean[i]/Q[i];
 171                                yres[i]  = ymean[i]-y;
 172
 173 // Systematic Error
 174 //
 175                                Int_t icx=segmentation.GetPadix(x);
 176                                Int_t icy=segmentation.GetPadiy(y);
 177                                xonpad[i]=segmentation.GetPadx(icx)-x;
 178                                yonpad[i]=segmentation.GetPady(icy)-y;
 179                            }
 180                        } // plane loop
 181                        xresid1->Fill(xres[0],(float) 1);
 182                        yresid1->Fill(yres[0],(float) 1);
 183                        npads1->Fill(npad[0],(float) 1);
 184                        if (npad[0] >=2 && Q[0] > 6 ) {
 185                            xresys1->Fill(xonpad[0],xres[0],(float) 1);
 186                            yresys1->Fill(yonpad[0],yres[0],(float) 1);
 187                        }
 188
 189                        xresid2->Fill(xres[1],(float) 1);
 190                        yresid2->Fill(yres[1],(float) 1);
 191                        npads2->Fill(npad[1],(float) 1);
 192                        if (npad[1] >=2 && Q[1] > 6) {
 193                            xresys2->Fill(xonpad[1],xres[1],(float) 1);
 194                            yresys2->Fill(yonpad[1],yres[1],(float) 1);
 195                        }
 196                    } // chamber 1/2
 197                } // hit loop
 198            } // if MUON
 199        } // track loop
 200    } // event loop
 201
 202 //Create a canvas, set the view range, show histograms
 203    TCanvas *c1 = new TCanvas("c1","Charge and Residuals",400,10,600,700);
 204    pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
 205    pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
 206    pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
 207    pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
 208    pad11->SetFillColor(11);
 209    pad12->SetFillColor(11);
 210    pad13->SetFillColor(11);
 211    pad14->SetFillColor(11);
 212    pad11->Draw();
 213    pad12->Draw();
 214    pad13->Draw();
 215    pad14->Draw();
 216
 217    pad11->cd();
 218    ccharge1->SetFillColor(42);
 219    ccharge1->SetXTitle("ADC units");
 220    ccharge1->Draw();
 221
 222    pad12->cd();
 223    pcharge1->SetFillColor(42);
 224    pcharge1->SetXTitle("ADC units");
 225    pcharge1->Draw();
 226
 227    pad13->cd();
 228    xresid1->SetFillColor(42);
 229    xresid1->Draw();
 230
 231    pad14->cd();
 232    yresid1->SetFillColor(42);
 233    yresid1->Draw();
 234
 235    TCanvas *c2 = new TCanvas("c2","Cluster Size",400,10,600,700);
 236    pad21 = new TPad("pad21"," ",0.01,0.51,0.49,0.99);
 237    pad22 = new TPad("pad22"," ",0.51,0.51,0.99,0.99);
 238    pad23 = new TPad("pad23"," ",0.01,0.01,0.49,0.49);
 239    pad24 = new TPad("pad24"," ",0.51,0.01,0.99,0.49);
 240    pad21->SetFillColor(11);
 241    pad22->SetFillColor(11);
 242    pad23->SetFillColor(11);
 243    pad24->SetFillColor(11);
 244    pad21->Draw();
 245    pad22->Draw();
 246    pad23->Draw();
 247    pad24->Draw();
 248
 249    pad21->cd();
 250    npads1->SetFillColor(42);
 251    npads1->SetXTitle("Cluster Size");
 252    npads1->Draw();
 253
 254    pad23->cd();
 255    xresys1->SetXTitle("x on pad");
 256    xresys1->SetYTitle("x-xcog");
 257    xresys1->Draw();
 258
 259    pad24->cd();
 260    yresys1->SetXTitle("y on pad");
 261    yresys1->SetYTitle("y-ycog");
 262    yresys1->Draw();
 263
 264    TCanvas *c3 = new TCanvas("c3","Charge and Residuals",400,10,600,700);
 265    pad31 = new TPad("pad31"," ",0.01,0.51,0.49,0.99);
 266    pad32 = new TPad("pad32"," ",0.51,0.51,0.99,0.99);
 267    pad33 = new TPad("pad33"," ",0.01,0.01,0.49,0.49);
 268    pad34 = new TPad("pad34"," ",0.51,0.01,0.99,0.49);
 269    pad31->SetFillColor(11);
 270    pad32->SetFillColor(11);
 271    pad33->SetFillColor(11);
 272    pad34->SetFillColor(11);
 273    pad31->Draw();
 274    pad32->Draw();
 275    pad33->Draw();
 276    pad34->Draw();
 277
 278    pad31->cd();
 279    ccharge2->SetFillColor(42);
 280    ccharge2->SetXTitle("ADC units");
 281    ccharge2->Draw();
 282
 283    pad32->cd();
 284    pcharge2->SetFillColor(42);
 285    pcharge2->SetXTitle("ADC units");
 286    pcharge2->Draw();
 287
 288    pad33->cd();
 289    xresid2->SetFillColor(42);
 290    xresid2->Draw();
 291
 292    pad34->cd();
 293    yresid2->SetFillColor(42);
 294    yresid2->Draw();
 295
 296    TCanvas *c4 = new TCanvas("c4","Cluster Size",400,10,600,700);
 297    pad41 = new TPad("pad41"," ",0.01,0.51,0.49,0.99);
 298    pad42 = new TPad("pad42"," ",0.51,0.51,0.99,0.99);
 299    pad43 = new TPad("pad43"," ",0.01,0.01,0.49,0.49);
 300    pad44 = new TPad("pad44"," ",0.51,0.01,0.99,0.49);
 301    pad41->SetFillColor(11);
 302    pad42->SetFillColor(11);
 303    pad43->SetFillColor(11);
 304    pad44->SetFillColor(11);
 305    pad41->Draw();
 306    pad42->Draw();
 307    pad43->Draw();
 308    pad44->Draw();
 309
 310    pad41->cd();
 311    npads2->SetFillColor(42);
 312    npads2->SetXTitle("Cluster Size");
 313    npads2->Draw();
 314
 315    pad43->cd();
 316    xresys2->SetXTitle("x on pad");
 317    xresys2->SetYTitle("x-xcog");
 318    xresys2->Draw();
 319
 320    pad44->cd();
 321    yresys2->SetXTitle("y on pad");
 322    yresys2->SetYTitle("y-ycog");
 323    yresys2->Draw();
 324
 325 }