MUON/MUONtestrawclust.C

   1 #include "iostream.h"
   2
   3 void MUONtestrawclust (Int_t evNumber1=0,Int_t evNumber2=0, Int_t ich1=0, Int_t ich2=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,-4,4);
  38    TH1F *hresy = new TH1F("hresy","Residuals",100,-.1,.1);
  39    TH1F *hresym = new TH1F("hresym","Residuals",100,-500,500);
  40    TH2F *hchi1 = new TH2F("hchi1","Chi2 vs Residuals",100,0,0.2,100,-500,500);
  41    TH2F *hchi2 = new TH2F("hchi2","Chi2 vs Residuals",100,0,20,100,-500,500);
  42 //
  43 //   Loop over events
  44 //
  45    Int_t Nh=0;
  46    Int_t Nh1=0;
  47    for (int nev=evNumber1; 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        AliMUONRawCluster  *mRaw;
  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        printf("Found %d entries in the tree (must be one per cathode per event! + 1empty)\n",nent);
  68        if (MUON) {
  69            for (Int_t ich=ich1;ich<=ich2;ich++) {
  70                TClonesArray *MUONrawclust  = MUON->RawClustAddress(ich);
  71                //          printf ("MUONrawclust %p \n",MUONrawclust);
  72
  73            for (Int_t icat=1; icat<2; icat++) {
  74                MUON->ResetRawClusters();
  75                nbytes += TR->GetEvent(icat);
  76                Int_t nrawcl = MUONrawclust->GetEntries();
  77                printf("Found %d raw clusters for cathode %d in chamber %d \n"
  78                       ,nrawcl,icat,ich+1);
  79                for (Int_t iraw=0; iraw < nrawcl; iraw++) {
  80                    mRaw = (AliMUONRawCluster*)MUONrawclust->UncheckedAt(iraw);
  81                    Int_t mult=mRaw->fMultiplicity[1];
  82                    Int_t itrack=mRaw->fTracks[1];
  83                    printf("\n mult1 mult2 %d %d chi2 %f itrack %d"
  84                           ,mRaw->fMultiplicity[0], mRaw->fMultiplicity[1], mRaw->fChi2[0], itrack);
  85                    h1->Fill(mult,float(1));
  86
  87                    Float_t xrec=mRaw->fX[1];
  88                    Float_t yrec=mRaw->fY[0];
  89                    Float_t R=TMath::Sqrt(xrec*xrec+yrec*yrec);
  90                    Int_t nres=0;
  91                    nbytes=0;
  92                    gAlice->ResetHits();
  93                    Int_t nbytes += TH->GetEvent(itrack);
  94
  95                        for(AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1);
  96                            mHit;
  97                            mHit=(AliMUONHit*)MUON->NextHit())
  98                        {
  99                            Int_t   nch   = mHit->fChamber;  // chamber number
 100                            Float_t x     = mHit->X();        // x-pos of hit
 101                            Float_t y     = mHit->Y();        // y-pos
 102                            if (nch==(ich+1)){
 103                                hresx->Fill(xrec-x,float(1));
 104                                hresy->Fill(yrec-y,float(1));
 105                                hchi1->Fill(mRaw->fChi2[0],(yrec-y)*1e4,float(1));
 106                                hchi2->Fill(mRaw->fChi2[0],(yrec-y)*1e4,float(1));
 107
 108                                if ((yrec-y)*1e4 <500 )
 109                                    hresym->Fill((yrec-y)*1e4,float(1));
 110                                if (mRaw->fChi2[0]>.3) {
 111                                    h22->Fill(mRaw->fX[1],mRaw->fY[0],float(1));
 112                                    hmult->Fill(mult,float(1));
 113                                }
 114                            } // chamber
 115                        } //hit
 116                } //iraw
 117            }  // icat
 118        } // ich
 119    }   // end if MUON
 120    }   // event loop
 121    TCanvas *c1 = new TCanvas("c1","Charge and Residuals",400,10,600,700);
 122    c1->Divide(2,2);
 123    c1->cd(1);
 124    hresx->SetFillColor(42);
 125    hresx->SetXTitle("xrec-x");
 126    hresx->Draw();
 127
 128    c1->cd(2);
 129    hresy->SetFillColor(42);
 130    hresy->SetXTitle("yrec-y");
 131    hresy->Draw();
 132
 133    c1->cd(3);
 134    h1->SetFillColor(42);
 135    h1->SetXTitle("multiplicity");
 136    h1->Draw();
 137
 138    c1->cd(4);
 139    hresym->SetFillColor(42);
 140    hresym->SetXTitle("yrec-y");
 141    hresym->Fit("gaus");
 142    hresym->Draw();
 143
 144    TCanvas *c2 = new TCanvas("c2","Charge and Residuals",400,10,600,700);
 145    c2->Divide(2,2);
 146
 147    c2->cd(1);
 148    h22->SetFillColor(42);
 149    h22->SetXTitle("x");
 150    h22->SetYTitle("y");
 151    h22->Draw();
 152
 153    c2->cd(2);
 154    hmult->SetFillColor(42);
 155    hmult->SetXTitle("multiplicity");
 156    hmult->Draw();
 157
 158    c2->cd(3);
 159    hchi1->SetFillColor(42);
 160    hchi1->SetXTitle("chi2");
 161    hchi1->Draw();
 162
 163    c2->cd(4);
 164    hchi2->SetFillColor(42);
 165    hchi2->SetXTitle("chi2");
 166    hchi2->Draw();
 167
 168 }
 169
 170
 171