MUON/MUONmultest.C

   1 #include "iostream.h"
   2
   3 void MUONmultest (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,-4,4);
  38    TH1F *hresy = new TH1F("hresy","Residuals",100,-.1,.1);
  39    TH1F *hresym = new TH1F("hresym","Residuals",100,-500,500);
  40 //
  41 //   Loop over events
  42 //
  43    Int_t Nh=0;
  44    Int_t Nh1=0;
  45    for (int nev=0; nev<= evNumber2; nev++) {
  46        Int_t nparticles = gAlice->GetEvent(nev);
  47        cout << "nev         " << nev <<endl;
  48        cout << "nparticles  " << nparticles <<endl;
  49        if (nev < evNumber1) continue;
  50        if (nparticles <= 0) return;
  51
  52        TTree *TH = gAlice->TreeH();
  53        Int_t ntracks = TH->GetEntries();
  54        cout<<"ntracks "<<ntracks<<endl;
  55
  56        Int_t nbytes = 0;
  57
  58        AliMUONRawCluster  *mRaw;
  59
  60 // Get pointers to Alice detectors and Digits containers
  61        AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
  62        TClonesArray *Particles = gAlice->Particles();
  63        TTree *TR = gAlice->TreeR();
  64        Int_t nent=TR->GetEntries();
  65        printf("Found %d entries in the tree (must be one per cathode per event! + 1empty)\n",nent);
  66        if (MUON) {
  67            for (Int_t ich=0;ich<1;ich++) {
  68                TClonesArray *MUONrawclust  = MUON->RawClustAddress(ich);
  69                TClonesArray *MUONdigits  = MUON->DigitsAddress(ich);
  70                for (Int_t icat=1; icat<2; icat++) {
  71                    MUON->ResetRawClusters();
  72                    nbytes += TR->GetEvent(icat);
  73                    Int_t nrawcl = MUONrawclust->GetEntries();
  74                    printf(
  75                        "Found %d raw-clusters for cathode %d in chamber %d \n"
  76                       ,nrawcl,icat,ich+1);
  77
  78
  79                    gAlice->ResetDigits();
  80
  81                    Int_t nent=(Int_t)gAlice->TreeD()->GetEntries();
  82                    gAlice->TreeD()->GetEvent(nent-2+icat-1);
  83                    Int_t ndigits = MUONdigits->GetEntriesFast();
  84                    printf(
  85                    "Found %d digits for cathode %d in chamber %d \n"
  86                        ,ndigits,icat,ich+1);
  87
  88
  89                    Int_t TotalMult =0;
  90
  91                    for (Int_t iraw=0; iraw < nrawcl; iraw++) {
  92                        mRaw = (AliMUONRawCluster*)MUONrawclust->UncheckedAt(iraw);
  93                        Int_t mult=mRaw->fMultiplicity;
  94                        h1->Fill(mult,float(1));
  95                        TotalMult+=mult;
  96                        Int_t itrack=mRaw->fTracks[0];
  97                        Float_t xrec=mRaw->fX;
  98                        Float_t yrec=mRaw->fY;
  99                        Float_t R=TMath::Sqrt(xrec*xrec+yrec*yrec);
 100                        if (R > 55.2) continue;
 101                        if (itrack ==1) continue;
 102                        Float_t res[2];
 103                        Int_t nres=0;
 104                        nbytes=0;
 105                        gAlice->ResetHits();
 106                        Int_t nbytes += TH->GetEvent(itrack);
 107
 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                            if (nch==(ich+1)){
 116                                hresx->Fill(xrec-x,float(1));
 117                                hresy->Fill(yrec-y,float(1));
 118                                if ((yrec-y)*1e4 <500 )
 119                                    hresym->Fill((yrec-y)*1e4,float(1));
 120                                if (TMath::Abs(yrec-y)>.02) {
 121                                    h22->Fill(mRaw->fX,mRaw->fY,float(1));
 122                                    hmult->Fill(mult,float(1));
 123                                }
 124                            } // chamber
 125                        } //hit
 126                    } //iraw
 127                    printf("Total Cluster Multiplicity %d \n" , TotalMult);
 128                }  // icat
 129            } // ich
 130        }   // end if MUON
 131    }   // event loop
 132    TCanvas *c1 = new TCanvas("c1","Charge and Residuals",400,10,600,700);
 133    pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
 134    pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
 135    pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
 136    pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
 137    pad11->SetFillColor(11);
 138    pad12->SetFillColor(11);
 139    pad13->SetFillColor(11);
 140    pad14->SetFillColor(11);
 141    pad11->Draw();
 142    pad12->Draw();
 143    pad13->Draw();
 144    pad14->Draw();
 145
 146    pad11->cd();
 147    hresx->SetFillColor(42);
 148    hresx->SetXTitle("xrec-x");
 149    hresx->Draw();
 150
 151    pad12->cd();
 152    hresy->SetFillColor(42);
 153    hresy->SetXTitle("yrec-y");
 154    hresy->Draw();
 155
 156    pad13->cd();
 157    h1->SetFillColor(42);
 158    h1->SetXTitle("multiplicity");
 159    h1->Draw();
 160
 161    pad14->cd();
 162    hresym->SetFillColor(42);
 163    hresym->SetXTitle("yrec-y");
 164    hresym->Fit("gaus");
 165    hresym->Draw();
 166    TCanvas *c2 = new TCanvas("c2","Charge and Residuals",400,10,600,700);
 167    pad21 = new TPad("pad21"," ",0.01,0.51,0.49,0.99);
 168    pad22 = new TPad("pad22"," ",0.51,0.51,0.99,0.99);
 169    pad23 = new TPad("pad23"," ",0.01,0.01,0.49,0.49);
 170    pad24 = new TPad("pad24"," ",0.51,0.01,0.99,0.49);
 171    pad21->SetFillColor(11);
 172    pad22->SetFillColor(11);
 173    pad23->SetFillColor(11);
 174    pad24->SetFillColor(11);
 175    pad21->Draw();
 176    pad22->Draw();
 177    pad23->Draw();
 178    pad24->Draw();
 179
 180    pad21->cd();
 181    h22->SetFillColor(42);
 182    h22->SetXTitle("x");
 183    h22->SetYTitle("y");
 184    h22->Draw();
 185
 186    pad22->cd();
 187    hmult->SetFillColor(42);
 188    hmult->SetXTitle("multiplicity");
 189    hmult->Draw();
 190
 191 }
 192
 193
 194