Added charge to the clusters tree
[u/mrichter/AliRoot.git] / ITS / UPGRADE / readHit.C
1 void readHit(int nev=-1,int evStart=0)
2 {
3   gROOT->SetStyle("Plain");
4
5   gSystem->Load("libITSUpgradeBase");
6   gSystem->Load("libITSUpgradeSim");
7
8   TH2F *xyGlob = new TH2F("xyGlob"," X - Y Global coordinates ",100,-50,50,100,-50,50);
9   xyGlob->SetXTitle("cm");
10   xyGlob->SetMarkerStyle(7);
11   TH1F *zGlob  = new TH1F("zGlob", " Z Global coordinates ",200, -100,100 );
12   zGlob->SetXTitle("cm");
13
14   Int_t nbins=100;
15   Double_t xmin=0;
16   Double_t xmax=1e-04;
17
18   gAlice=NULL;
19   AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
20   runLoader->LoadgAlice();
21
22   gAlice = runLoader->GetAliRun();
23
24   runLoader->LoadHeader();
25   runLoader->LoadKinematics();
26   runLoader->LoadSDigits();
27   runLoader->LoadHits();
28
29   AliGeomManager::LoadGeometry("geometry.root");
30   AliITSUGeomTGeo* gm = new AliITSUGeomTGeo(kTRUE);
31
32   AliLoader *dl = runLoader->GetDetectorLoader("ITS");
33
34
35
36   Int_t nLayers = gm->GetNLayers();
37
38   TH1D **hDeLoss = new TH1D*[nLayers];
39   for(Int_t i=0; i< nLayers; i++ ) {
40     hDeLoss[i] = new TH1D(Form("hDeLossl%i",i),Form("E loss distribution [ Layer %i] ",i),nbins,xmin,xmax);
41     hDeLoss[i]->SetXTitle("GeV");
42   }
43
44   int nevTot = (Int_t)runLoader->GetNumberOfEvents();
45   printf("N Events : %i \n",nevTot);
46   evStart = evStart<nevTot ? evStart : nevTot-1;
47   if (evStart<0) evStart = 0;
48   //
49   int lastEv = nev<0 ? nevTot : evStart+nev;
50   if (lastEv > nevTot) lastEv = nevTot;
51   //
52
53   //HIT INIT
54
55   TTree *hitTree = 0x0;
56   TClonesArray *hitList=new TClonesArray("AliITSUHit");
57
58   for (Int_t iEvent = evStart; iEvent < lastEv; iEvent++) {
59
60     printf("\nEvent\t%d\n",iEvent);
61  
62     runLoader->GetEvent(iEvent);
63     
64     hitTree=dl->TreeH();
65     hitTree->SetBranchAddress("ITS",&hitList);
66     for(Int_t iEnt=0;iEnt<hitTree->GetEntries();iEnt++){//entries loop degli hits
67       hitTree->GetEntry(iEnt);
68       for(Int_t iHit=0; iHit<hitList->GetEntries();iHit++){
69               
70         AliITSUHit *pHit = (AliITSUHit*)hitList->At(iHit);
71         int id = pHit->GetModule();
72         int lr = gm->GetLayer(id);
73         int ld = gm->GetLadder(id);
74         //
75         //      if(pHit->GetParticle()->IsPrimary()){
76         Double_t xg,yg,zg=0.;
77         Double_t xg0,yg0,zg0=0.,tg0;
78         pHit->GetPositionG(xg,yg,zg);
79         pHit->GetPositionG0(xg0,yg0,zg0,tg0);
80         xyGlob->Fill(xg,yg);
81         zGlob->Fill(zg);
82         printf("Module %5d | Lr:%2d Ladder: %3d, X:[%+.5e:%+.5e] Y:[%+.5e:%+.5e] Z:[%+.5e %+.5e] DE: %.5e TrackID: %d\n",id,lr,ld,
83                xg0,xg,yg0,yg,zg0,zg,pHit->GetIonization(),pHit->GetTrack());
84         hDeLoss[lr]->Fill(pHit->GetIonization());
85         //      } // is primary
86       }//loop hit 
87     }//entryloopHitList
88         
89   }//event loop
90
91   TCanvas *xyCanv =  new TCanvas("xyCanv","Hit X-Y positions",500,500);
92   xyCanv->cd();
93   xyGlob->Draw();
94   TCanvas *zCanv =  new TCanvas("zCanv","Hit Z positions",500,500);
95   zCanv->cd();
96   zGlob->Draw();
97
98   TCanvas *c = new TCanvas("c","E loss  distribution per layer",1000,800);
99   c->Divide(nLayers&0x1 ?  (nLayers/2+1) : nLayers/2,2);
100   for(Int_t ip =1; ip<=nLayers; ip++){
101     c->cd(ip);
102     hDeLoss[ip-1]->Draw();
103   }
104
105 }