MUON/MUONtestsim.C

   1 void MUONtestsim(Int_t ncham, 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
  13     if (gClassTable->GetID("AliRun") < 0) {
  14         gROOT->LoadMacro("loadlibs.C");
  15         loadlibs();
  16     }
  17
  18 // Connect the Root Galice file containing Geometry, Kine and Hits
  19
  20     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  21
  22
  23     if (!file) {
  24         printf("\n Creating galice.root \n");
  25         file = new TFile("galice.root");
  26     } else {
  27         printf("\n galice.root found in file list");
  28     }
  29     file->ls();
  30
  31 // Get AliRun object from file or create it if not on file
  32     if (!gAlice) {
  33         gAlice = (AliRun*)(file->Get("gAlice"));
  34         if (gAlice) printf("AliRun object found on file\n");
  35         if (!gAlice) {
  36             printf("\n create new gAlice object");
  37             gAlice = new AliRun("gAlice","Alice test program");
  38         }
  39     }
  40
  41 //  Create some histograms
  42
  43
  44    TH1F *ccharge1 = new TH1F("ccharge1","Cluster Charge distribution"
  45                            ,100,0.,5000.);
  46    TH1F *pcharge1 = new TH1F("pcharge1","Pad      Charge distribution"
  47                            ,100,0.,200.);
  48    TH1F *xresid1  = new TH1F("xresid1","x-Residuals"
  49                            ,100,-10.,10);
  50    TH1F *yresid1  = new TH1F("yresid1","y-Residuals"
  51                            ,100,-0.2,0.2);
  52    TH1F *npads1   = new TH1F("npads1" ,"Pads per Hit"
  53                            ,20,-0.5,19.5);
  54    TH2F *xresys1  = new TH2F("xresys1","x-Residuals systematics"
  55                            ,50,-2.0,2.0,100,-1.0,1.0);
  56    TH2F *yresys1  = new TH2F("yresys1","y-Residuals systematics"
  57                            ,50,-0.4,0.4,100,-0.1,0.1);
  58
  59    TH1F *ccharge2 = new TH1F("ccharge2","Cluster Charge distribution"
  60                            ,100,0.,5000.);
  61    TH1F *pcharge2 = new TH1F("pcharge2","Pad      Charge distribution"
  62                            ,100,0.,200.);
  63    TH1F *xresid2  = new TH1F("xresid2","x-Residuals"
  64                            ,100,-1,1);
  65    TH1F *yresid2  = new TH1F("yresid2","y-Residuals"
  66                            ,100,-10, 10.);
  67    TH1F *npads2   = new TH1F("npads2" ,"Pads per Hit"
  68                            ,20,-0.5,19.5);
  69    TH2F *xresys2  = new TH2F("xresys2","x-Residuals systematics"
  70                            ,50,-2.0,2.0,100,-1.0,1.0);
  71    TH2F *yresys2  = new TH2F("yresys2","y-Residuals systematics"
  72                            ,50,-0.4,0.4,100,-0.1,0.1);
  73
  74
  75    AliMUONChamber*  iChamber;
  76 //
  77 //   Loop over events
  78 //
  79    Int_t Nh=0;
  80    Int_t Nh1=0;
  81    for (Int_t nev=0; nev<= evNumber2; nev++) {
  82        cout << "nev         " << nev <<endl;
  83        Int_t nparticles = gAlice->GetEvent(nev);
  84        cout << "nparticles  " << nparticles <<endl;
  85        if (nev < evNumber1) continue;
  86        if (nparticles <= 0) return;
  87
  88        AliMUON *MUON  = (AliMUON*) gAlice->GetDetector("MUON");
  89            printf("\n track %d %d \n ", nev, MUON);
  90
  91        TTree *TH = gAlice->TreeH();
  92        Int_t ntracks = TH->GetEntries();
  93        Int_t Nc=0;
  94        Int_t   npad[2];
  95        Float_t Q[2], xmean[2],ymean[2],xres[2],yres[2], xonpad[2], yonpad[2];
  96 //
  97 //   Loop over events
  98 //
  99        for (Int_t track=0; track<ntracks;track++) {
 100
 101            gAlice->ResetHits();
 102            Int_t nbytes += TH->GetEvent(track);
 103            if (MUON)  {
 104                for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
 105                    mHit;
 106                    mHit=(AliMUONHit*)MUON->NextHit())
 107                {
 108                    Int_t   nch   = mHit->fChamber;  // chamber number
 109                    Float_t x     = mHit->X();        // x-pos of hit
 110                    Float_t y     = mHit->Y();        // y-pos
 111                    Float_t z     = mHit->Z();
 112
 113 //
 114 //
 115                    iChamber = & MUON->Chamber(nch-1);
 116                    response=iChamber->ResponseModel();
 117 //
 118 //
 119                    if (nch == ncham) {
 120                        for (Int_t i=0; i<2; i++) {
 121                            xmean[i]=0;
 122                            ymean[i]=0;
 123                            Q[i]=0;
 124                            npad[i]=0;
 125                        }
 126
 127                        for (AliMUONPadHit* mPad=
 128                                 (AliMUONPadHit*)MUON->FirstPad(mHit,MUON->PadHits());
 129                             mPad;
 130                             mPad=(AliMUONPadHit*)MUON->NextPad(MUON->PadHits()))
 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 secpad   = mPad->fRSec;      // r-pos of pad
 140 //
 141 //
 142                            segmentation=iChamber->SegmentationModel(nseg);
 143                            Int_t ind=nseg-1;
 144                            Float_t xg,yg,zg;
 145                            segmentation->GetPadC(ipx, ipy, xg, yg, zg);
 146                            Int_t sec = segmentation->Sector(ipx, ipy);
 147                            Float_t dpx=segmentation->Dpx(sec);
 148                            Float_t dpy=segmentation->Dpy(sec);
 149
 150                            if (nseg == 1) {
 151                                pcharge1->Fill(iqpad,(float) 1);
 152                                ccharge1->Fill(qtot ,(float) 1);
 153                            }  else {
 154                                pcharge2->Fill(iqpad,(float) 1);
 155                                ccharge2->Fill(qtot ,(float) 1);
 156                            }
 157 // Calculate centre of gravity
 158 //
 159                            if (iqpad == 0 && ind==0) {
 160                                printf("\n attention iqpad 0");
 161                            }
 162
 163                            if (iqpad >  0) {
 164                                Float_t xc,yc,zc;
 165                                npad[ind]++;
 166                                segmentation->GetPadC(ipx,ipy,xc,yc,zc);
 167                                xmean[ind]+=Float_t(iqpad*xc);
 168                                ymean[ind]+=Float_t(iqpad*yc);
 169                                Q[ind]+=iqpad;
 170                            }
 171
 172                        } //pad hit loop
 173                        for (Int_t i=0; i<2; i++) {
 174                            segmentation = iChamber->SegmentationModel(i+1);
 175                            if (Q[i] >0) {
 176                                xmean[i] = xmean[i]/Q[i];
 177                                xres[i]  = xmean[i]-x;
 178                                ymean[i] = ymean[i]/Q[i];
 179                                yres[i]  = ymean[i]-y;
 180 // Systematic Error
 181 //
 182                                Int_t icx, icy;
 183                                Float_t zonpad;
 184
 185                                segmentation->
 186                                    GetPadI(x,y,z,icx,icy);
 187                                segmentation->
 188                                    GetPadC(icx,icy,xonpad[i],yonpad[i],zonpad);
 189                                xonpad[i]-=x;
 190                                yonpad[i]-=y;
 191                            } // charge not 0
 192                        } // plane loop
 193                        xresid1->Fill(xres[0],(float) 1);
 194                        yresid1->Fill(yres[0],(float) 1);
 195                        npads1->Fill(npad[0],(float) 1);
 196                        if (npad[0] >=2 && Q[0] > 6 ) {
 197                            xresys1->Fill(xonpad[0],xres[0],(float) 1);
 198                            yresys1->Fill(yonpad[0],yres[0],(float) 1);
 199                        }
 200
 201                        xresid2->Fill(xres[1],(float) 1);
 202                        yresid2->Fill(yres[1],(float) 1);
 203                        npads2->Fill(npad[1],(float) 1);
 204                        if (npad[1] >=2 && Q[1] > 6) {
 205                            xresys2->Fill(xonpad[1],xres[1],(float) 1);
 206                            yresys2->Fill(yonpad[1],yres[1],(float) 1);
 207                        }
 208                    } // chamber 1
 209                } // hit loop
 210            } // if MUON
 211        } // track loop
 212    } // event loop
 213
 214 //Create a canvas, set the view range, show histograms
 215    TCanvas *c1 = new TCanvas("c1","Charge and Residuals (1)",400,10,600,700);
 216    pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
 217    pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
 218    pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
 219    pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
 220    pad11->SetFillColor(11);
 221    pad12->SetFillColor(11);
 222    pad13->SetFillColor(11);
 223    pad14->SetFillColor(11);
 224    pad11->Draw();
 225    pad12->Draw();
 226    pad13->Draw();
 227    pad14->Draw();
 228
 229    pad11->cd();
 230    ccharge1->SetFillColor(42);
 231    ccharge1->SetXTitle("ADC units");
 232    ccharge1->Draw();
 233
 234    pad12->cd();
 235    pcharge1->SetFillColor(42);
 236    pcharge1->SetXTitle("ADC units");
 237    pcharge1->Draw();
 238
 239    pad13->cd();
 240    xresid1->SetFillColor(42);
 241    xresid1->Draw();
 242
 243    pad14->cd();
 244    yresid1->SetFillColor(42);
 245    yresid1->Draw();
 246
 247    TCanvas *c2 = new TCanvas("c2","Cluster Size (1)",400,10,600,700);
 248    pad21 = new TPad("pad21"," ",0.01,0.51,0.49,0.99);
 249    pad22 = new TPad("pad22"," ",0.51,0.51,0.99,0.99);
 250    pad23 = new TPad("pad23"," ",0.01,0.01,0.49,0.49);
 251    pad24 = new TPad("pad24"," ",0.51,0.01,0.99,0.49);
 252    pad21->SetFillColor(11);
 253    pad22->SetFillColor(11);
 254    pad23->SetFillColor(11);
 255    pad24->SetFillColor(11);
 256    pad21->Draw();
 257    pad22->Draw();
 258    pad23->Draw();
 259    pad24->Draw();
 260
 261    pad21->cd();
 262    npads1->SetFillColor(42);
 263    npads1->SetXTitle("Cluster Size");
 264    npads1->Draw();
 265
 266    pad23->cd();
 267    xresys1->SetXTitle("x on pad");
 268    xresys1->SetYTitle("x-xcog");
 269    xresys1->Draw();
 270
 271    pad24->cd();
 272    yresys1->SetXTitle("y on pad");
 273    yresys1->SetYTitle("y-ycog");
 274    yresys1->Draw();
 275
 276    TCanvas *c3 = new TCanvas("c3","Charge and Residuals (2)",400,10,600,700);
 277    pad31 = new TPad("pad31"," ",0.01,0.51,0.49,0.99);
 278    pad32 = new TPad("pad32"," ",0.51,0.51,0.99,0.99);
 279    pad33 = new TPad("pad33"," ",0.01,0.01,0.49,0.49);
 280    pad34 = new TPad("pad34"," ",0.51,0.01,0.99,0.49);
 281    pad31->SetFillColor(11);
 282    pad32->SetFillColor(11);
 283    pad33->SetFillColor(11);
 284    pad34->SetFillColor(11);
 285    pad31->Draw();
 286    pad32->Draw();
 287    pad33->Draw();
 288    pad34->Draw();
 289
 290    pad31->cd();
 291    ccharge2->SetFillColor(42);
 292    ccharge2->SetXTitle("ADC units");
 293    ccharge2->Draw();
 294
 295    pad32->cd();
 296    pcharge2->SetFillColor(42);
 297    pcharge2->SetXTitle("ADC units");
 298    pcharge2->Draw();
 299
 300    pad33->cd();
 301    xresid2->SetFillColor(42);
 302    xresid2->Draw();
 303
 304    pad34->cd();
 305    yresid2->SetFillColor(42);
 306    yresid2->Draw();
 307
 308    TCanvas *c4 = new TCanvas("c4","Cluster Size (2)",400,10,600,700);
 309    pad41 = new TPad("pad41"," ",0.01,0.51,0.49,0.99);
 310    pad42 = new TPad("pad42"," ",0.51,0.51,0.99,0.99);
 311    pad43 = new TPad("pad43"," ",0.01,0.01,0.49,0.49);
 312    pad44 = new TPad("pad44"," ",0.51,0.01,0.99,0.49);
 313    pad41->SetFillColor(11);
 314    pad42->SetFillColor(11);
 315    pad43->SetFillColor(11);
 316    pad44->SetFillColor(11);
 317    pad41->Draw();
 318    pad42->Draw();
 319    pad43->Draw();
 320    pad44->Draw();
 321
 322    pad41->cd();
 323    npads2->SetFillColor(42);
 324    npads2->SetXTitle("Cluster Size");
 325    npads2->Draw();
 326
 327    pad43->cd();
 328    xresys2->SetXTitle("x on pad");
 329    xresys2->SetYTitle("x-xcog");
 330    xresys2->Draw();
 331
 332    pad44->cd();
 333    yresys2->SetXTitle("y on pad");
 334    yresys2->SetYTitle("y-ycog");
 335    yresys2->Draw();
 336
 337 }