MUON/MUONhitrectest.C

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