Modifications needed by the HBT analysis (P.Skowronski)
[u/mrichter/AliRoot.git] / MUON / MUONtest.C
1 void MUONtestabso (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         gSystem->Load("$ALITOP/cern.so/lib/libpdfDUMMY.so");
14         gSystem->Load("$ALITOP/cern.so/lib/libPythia.so");
15         gSystem->Load("$ROOTSYS/lib/libEG.so");       
16         gSystem->Load("$ROOTSYS/lib/libEGPythia.so");       
17         gSystem->Load("libGeant3Dummy.so");    //a dummy version of Geant3
18         gSystem->Load("PHOS/libPHOSdummy.so"); //the standard Alice classes 
19         gSystem->Load("libgalice.so");         // the standard Alice classes 
20     }
21 // Connect the Root Galice file containing Geometry, Kine and Hits
22
23     TFile *file = (TFile*)gROOT->GetListOfFiles()->FindObject("galice.root");
24
25     
26     if (!file) {
27         printf("\n Creating galice.root \n");
28         file = new TFile("galice.root");
29     } else {
30         printf("\n galice.root found in file list");
31     }
32     file->ls();
33     
34 // Get AliRun object from file or create it if not on file
35     if (!gAlice) {
36         gAlice = (AliRun*)(file->Get("gAlice"));
37         if (gAlice) printf("AliRun object found on file\n");
38         if (!gAlice) {
39             printf("\n create new gAlice object");
40             gAlice = new AliRun("gAlice","Alice test program");
41         }
42     }
43     
44 //  Create some histograms
45
46
47    TH1F *zv = new TH1F("zv","z vertex"
48                            ,100,468.,503.);
49    TH1F *xv = new TH1F("xv","x vertex"
50                            ,100,0.,50.);
51    TH1F *yv  = new TH1F("yv","y vertex"
52                            ,100,0.,50.);
53
54    TH1F *ip  = new TH1F("ip","ipart"
55                            ,100,0.,50.);
56
57    AliMUONchamber*  iChamber;
58 //
59 //   Loop over events 
60 //
61    Int_t Nh=0;
62    Int_t Nh1=0;
63    for (Int_t nev=0; nev<= evNumber2; nev++) {
64        cout << "nev         " << nev <<endl;
65        Int_t nparticles = gAlice->GetEvent(nev);
66        cout << "nparticles  " << nparticles <<endl;
67        if (nev < evNumber1) continue;
68        if (nparticles <= 0) return;
69
70        AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
71            printf("\n nev %d \n ", nev);
72
73        TTree *TH = gAlice->TreeH();
74        Int_t ntracks = TH->GetEntries();
75 //
76 //   Loop over tracks
77 //
78        for (Int_t track=0; track<ntracks;track++) {
79            
80            gAlice->ResetHits();
81            Int_t nbytes += TH->GetEvent(track);
82            if (MUON)  {
83                for(AliMUONhit* mHit=(AliMUONhit*)MUON->FirstHit(-1); 
84                    mHit;
85                    mHit=(AliMUONhit*)MUON->NextHit()) 
86                {
87                    Int_t   nch   = mHit->fChamber;  // chamber number
88                    Float_t x     = mHit->fX;        // x-pos of hit
89                    Float_t y     = mHit->fY;        // y-pos
90                    Float_t z     = mHit->fZ;        // y-pos
91                
92                 if (nch > 1) continue;
93
94                 Int_t ipart = mHit->fParticle;
95                 TClonesArray *fPartArray = gAlice->Particles();
96                 TParticle *Part;
97                 Int_t ftrack = mHit->fTrack;
98                 Int_t id = ((TParticle*) fPartArray->UncheckedAt(ftrack))->GetPdgCode();
99                   ip->Fill((float)ipart,(float) 1);
100
101                 if (ipart != 3) continue;
102                   Float_t xv = Part->Vx();      //  vertex  
103                   Float_t yv = Part->Vy();      //  vertex
104                   Float_t zv  = Part->Vz();      // z vertex 
105                   xv->Fill(xv,(float) 1);
106                   yv->Fill(yv,(float) 1);
107                   zv->Fill(zv,(float) 1);
108                }
109
110            }
111        }
112    }
113
114    
115 //Create a canvas, set the view range, show histograms
116    TCanvas *c1 = new TCanvas("c1","Charge and Residuals",400,10,600,700);
117    pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
118    pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
119    pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
120    pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
121    pad11->SetFillColor(11);
122    pad12->SetFillColor(11);
123    pad13->SetFillColor(11);
124    pad14->SetFillColor(11);
125    pad11->Draw();
126    pad12->Draw();
127    pad13->Draw();
128    pad14->Draw();
129    
130    pad11->cd();
131    ip->SetFillColor(42);
132    ip->SetXTitle("ipart");
133    ip->Draw();
134
135    pad12->cd();
136    xv->SetFillColor(42);
137    xv->SetXTitle("xvert");
138    xv->Draw();
139
140    pad13->cd();
141    yv->SetFillColor(42);
142    yv->SetXTitle("yvert");
143    yv->Draw();
144
145    pad14->cd();
146    zv->SetFillColor(42);
147    zv->SetXTitle("zvert");
148    zv->Draw();
149    
150 }