MUON/MUONhitrectest.C

   1 void MUONhitrectest (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       gROOT->LoadMacro("loadlibs.C");
  14       loadlibs();
  15    }
  16
  17 // Connect the Root Galice file containing Geometry, Kine and Hits
  18
  19    TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
  20    if (!file) file = new TFile("galice.root");
  21
  22 // Get AliRun object from file or create it if not on file
  23
  24    if (!gAlice) {
  25       gAlice = (AliRun*)file->Get("gAlice");
  26       if (gAlice) printf("AliRun object found on file\n");
  27       if (!gAlice) gAlice = new AliRun("gAlice","Alice test program");
  28    }
  29 //  Create some histograms
  30
  31    TH2F *h21 = new TH2F("h21","Hits",100,-100,100,100,-100,100);
  32    TH2F *h22 = new TH2F("h22","CoG ",100,-100,100,100,-100,100);
  33    TH1F *h1 = new TH1F("h1","Multiplicity",30,-0.5,29.5);
  34    TH1F *hmult = new TH1F("hmult","Multiplicity",30,-0.5,29.5);
  35    TH1F *hresx = new TH1F("hresx","Residuals",100,-1,1);
  36    TH1F *hresy = new TH1F("hresy","Residuals",100,-10.,10);
  37    TH1F *hresym = new TH1F("hresym","Residuals",100,-500,500);
  38    TH1F *hpos = new TH1F("hnpos","Possibilities",10,-0.5,9.5);
  39    TH2F *hchi1 = new TH2F("hchi1","Chi2 vs Residuals",100,0,0.2,100,-500,500);
  40    TH2F *hchi2 = new TH2F("hchi2","Chi2 vs Residuals",100,0,20,100,-500,500);
  41 //
  42 //   Loop over events
  43 //
  44    Int_t Nh=0;
  45    Int_t Nh1=0;
  46
  47    for (int nev=0; nev<= evNumber2; nev++) {
  48        Int_t nparticles = gAlice->GetEvent(nev);
  49        cout << "nev         " << nev <<endl;
  50        cout << "nparticles  " << nparticles <<endl;
  51        if (nev < evNumber1) continue;
  52        if (nparticles <= 0) return;
  53
  54        TTree *TH = gAlice->TreeH();
  55        Int_t ntracks = TH->GetEntries();
  56        cout<<"ntracks "<<ntracks<<endl;
  57
  58        Int_t nbytes = 0;
  59
  60
  61
  62 // Get pointers to Alice detectors and Digits containers
  63        AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
  64        TClonesArray *Particles = gAlice->Particles();
  65        TTree *TR = gAlice->TreeR();
  66        Int_t nent=TR->GetEntries();
  67
  68        if (MUON) {
  69            Float_t xc[2000], yc[2000], itrc[2000];
  70            TClonesArray *MUONrawclust  = MUON->RawClustAddress(0);
  71            nbytes += TR->GetEvent(0);
  72            Int_t nrawcl = MUONrawclust->GetEntries();
  73            AliMUONRawCluster  *mRaw;
  74            for (Int_t iraw=0; iraw < nrawcl; iraw++) {
  75                mRaw = (AliMUONRawCluster*)MUONrawclust->UncheckedAt(iraw);
  76                xc[iraw]=mRaw->fX[1];
  77                yc[iraw]=mRaw->fY[0];
  78                itrc[iraw]=mRaw->fTracks[1];
  79            } // cluster
  80
  81            gAlice->ResetHits();
  82            for (Int_t itrack=0; itrack<ntracks ;itrack++)
  83            {
  84                Int_t nbytes += TH->GetEvent(itrack);
  85                Int_t nhit=MUON->Hits()->GetEntriesFast();
  86                for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(itrack);
  87                    mHit;
  88                    mHit=(AliMUONHit*)MUON->NextHit())
  89                {
  90                    Int_t   nch   = mHit->Chamber();  // chamber number
  91                    Float_t x     = mHit->X();        // x-pos of hit
  92                    Float_t y     = mHit->Y();        // y-pos
  93                    Int_t   ip    = mHit->Particle();
  94                    if (ip != kMuonPlus && ip != kMuonMinus) continue;
  95                    if (nch != 1) continue;
  96                    Float_t dmin=1.e8;
  97                    Int_t imin=-1;
  98 //
  99 // Loop over clusters
 100                    Int_t npos=0;
 101
 102                    for (Int_t i=0; i<nrawcl; i++) {
 103                        Float_t dx = xc[i]-x;
 104                        Float_t dy = yc[i]-y;
 105                        Float_t dr = TMath::Sqrt(dx*dx+dy*dy);
 106                        if (dr<dmin) {
 107                            imin=i;
 108                            dmin=dr;
 109                        }
 110                        if (TMath::Abs(dy) < 0.05 && TMath::Abs(dx) < 0.2 ) npos++;
 111                    }
 112                    mRaw = (AliMUONRawCluster*)MUONrawclust->UncheckedAt(imin);
 113                    Int_t   track = mRaw->fTracks[1];
 114                    Float_t xrec  = mRaw->fX[1];
 115                    Float_t yrec  = mRaw->fY[0];
 116                    Float_t x_res = xrec-x;
 117                    Float_t y_res = yrec-y;
 118                    if (y_res > 0.05 ) printf("\n track %d %f\n", itrack, y_res);
 119
 120                    hresx->Fill(x_res,1.);
 121                    hresy->Fill(y_res,1.);
 122                    hpos->Fill(Float_t(npos), 1.);
 123                    hresym->Fill(y_res*1.e4,1.);
 124                    if (npos == 0) {
 125                        printf("No cluster %d %f %f %d\n",
 126                               itrack, x, y, nhit);
 127                        h21->Fill(x,y,1.);
 128                    }
 129                }  // hits
 130            } // tracks
 131        }   // end if MUON
 132    }   // event loop
 133    TCanvas *c1 = new TCanvas("c1","Charge and Residuals",400,10,600,700);
 134    c1->Divide(2,2);
 135    c1->cd(1);
 136    hresx->SetFillColor(42);
 137    hresx->SetXTitle("xrec-x");
 138    hresx->Draw();
 139
 140    c1->cd(2);
 141    hresy->SetFillColor(42);
 142    hresy->SetXTitle("yrec-y");
 143    hresy->Draw();
 144
 145    c1->cd(3);
 146    hpos->SetFillColor(42);
 147    hpos->SetXTitle("Possibilities");
 148    hpos->Draw();
 149
 150    c1->cd(4);
 151    hresym->SetFillColor(42);
 152    hresym->SetXTitle("yrec-y");
 153    hresym->Fit("gaus");
 154    hresym->Draw();
 155
 156    TCanvas *c2 = new TCanvas("c2","Charge and Residuals",400,10,600,700);
 157    c2->Divide(2,2);
 158
 159    c2->cd(1);
 160    h21->SetFillColor(42);
 161    h21->SetXTitle("x");
 162    h21->SetYTitle("y");
 163    h21->Draw();
 164
 165    c2->cd(2);
 166    hmult->SetFillColor(42);
 167    hmult->SetXTitle("multiplicity");
 168    hmult->Draw();
 169
 170    c2->cd(3);
 171    hchi1->SetFillColor(42);
 172    hchi1->SetXTitle("chi2");
 173    hchi1->Draw();
 174
 175    c2->cd(4);
 176    hchi2->SetFillColor(42);
 177    hchi2->SetXTitle("chi2");
 178    hchi2->Draw();
 179
 180 }
 181
 182
 183