Added to TrackerGlo possibility to fill control histos with residuals,pulls,chi2
[u/mrichter/AliRoot.git] / ITS / UPGRADE / readDigits.C
1 void readDigits(){
2
3   gSystem->Load("libITSUpgradeBase");
4   gSystem->Load("libITSUpgradeSim");
5   gROOT->SetStyle("Plain");
6   const Int_t kMaxROCycleAccept=126;
7   gAlice=NULL;
8   AliRunLoader* runLoader = AliRunLoader::Open("galice.root");
9   runLoader->LoadgAlice();
10
11   gAlice = runLoader->GetAliRun();
12
13   runLoader->LoadHeader();
14   runLoader->LoadKinematics();
15   runLoader->LoadSDigits();
16   runLoader->LoadDigits();
17
18   AliGeomManager::LoadGeometry("geometry.root");
19   AliITSUGeomTGeo* gm = new AliITSUGeomTGeo(kTRUE,kTRUE);
20   //
21   Int_t nLayers = gm->GetNLayers();
22   Int_t nModules = gm->GetNModules();
23   Int_t nbins=100;
24   Int_t xmin=0;
25   Int_t xmax=50000;//00*1e-09;
26
27
28   TH1D **hNelSDig = new TH1D*[nLayers], **hNelDig = new TH1D*[nLayers];
29   //
30   for(Int_t i=0; i< nLayers; i++ ) {
31     hNelSDig[i] = new TH1D(Form("hNelSDig%i",i),Form("electron distribution in SDigits [ Layer %i] ",i),nbins,xmin,xmax);
32     hNelSDig[i]->SetXTitle("N electrons");
33     hNelDig[i] = new TH1D(Form("hNelDig%i",i),Form("electron distribution in Digits [ Layer %i] ",i),nbins,xmin,xmax);
34     hNelDig[i]->SetXTitle("N electrons");
35   }
36   
37   AliLoader *dl = runLoader->GetDetectorLoader("ITS");
38
39
40   //SDIGITS INIT
41   TTree * sDigTree = 0x0;
42   TClonesArray *sDigArr= new TClonesArray("AliITSUSDigit");
43
44   //DIGITS INIT
45   TTree * digTree = 0x0;
46   TClonesArray *digArr= new TClonesArray("AliITSUDigitPix");
47
48   printf("N Events : %i \n",(Int_t)runLoader->GetNumberOfEvents());
49
50   for (Int_t iEvent = 0; iEvent < runLoader->GetNumberOfEvents(); iEvent++) {
51     printf("\n Event %i \n",iEvent);
52     runLoader->GetEvent(iEvent);
53     AliStack *stack = runLoader->Stack();
54     sDigTree=dl->TreeS();
55     digTree=dl->TreeD();
56     //
57     sDigTree->SetBranchAddress("ITS",&sDigArr);
58     digTree->SetBranchAddress("ITSDigitsPix",&digArr);
59
60     for (int imod=0;imod<nModules;imod++) {
61       sDigTree->GetEntry(imod);
62       digTree->GetEntry(imod);      
63       int detType = gm->GetModuleDetTypeID(imod);
64       AliITSUSegmentationPix* segm = (AliITSUSegmentationPix*)gm->GetSegmentationByID(detType);
65       int lay,lad,det;
66       int nsdig = sDigArr->GetEntries();
67       int ndig  = digArr->GetEntries();
68       if (nsdig<1) continue;
69       gm->GetModuleId(imod, lay,lad,det);
70       printf("\nModule %3d: (det %2d in ladder %2d of Layer %d) |NSDigits: %4d NDigits: %4d\n",imod,det,lad,lay,nsdig,ndig);
71       //
72       for (int isdig=0;isdig<nsdig;isdig++) {
73         AliITSUSDigit *pSdig = (AliITSUSDigit*)sDigArr->At(isdig);
74         int sdinfo = pSdig->GetUniqueID();
75         UInt_t row,col;
76         Int_t cycle;
77         AliITSUSensMap::GetCell(sdinfo,segm->Npz(),segm->Npx(),kMaxROCycleAccept,col,row,cycle);
78         printf("#%3d Sdigit col:%4d/row:%4d/cycle:%d generated by track %5d (%s)\t",isdig, col,row,cycle,
79                pSdig->GetTrack(0),stack->Particle(pSdig->GetTrack(0))->GetName());
80         pSdig->Print();
81         hNelSDig[lay]->Fill(pSdig->GetSignal()); 
82       }
83       //
84       for (int idig=0;idig<ndig;idig++) {
85         AliITSUDigitPix *pDig = (AliITSUDigitPix*)digArr->At(idig);
86         printf("#%3d digit, col:%4d/row:%4d ROCycle:%d signal: %.2e,  generated by tracks ",idig,pDig->GetCoord1(),pDig->GetCoord2(),
87                pDig->GetROCycle(),pDig->GetSignalPix()); 
88         for (int itr=0;itr<pDig->GetNTracks();itr++) if (pDig->GetTrack(itr)>=0) printf(" %5d",pDig->GetTrack(itr)); printf("\n");
89         //
90         hNelDig[lay]->Fill(pDig->GetSignalPix()); 
91       }
92       //
93     }
94
95   }//event loop
96
97   TCanvas *cSd = new TCanvas("cSd","Summable Digits",1000,800);
98   cSd->Divide(nLayers&0x1 ?  (nLayers/2+1) : nLayers/2,2); 
99   TCanvas *cD = new TCanvas("cD","Digits",1000,800);
100   cD->Divide(nLayers&0x1 ?  (nLayers/2+1) : nLayers/2,2); 
101
102   for(Int_t ip =1; ip<=nLayers; ip++){
103     cSd->cd(ip);
104     hNelSDig[ip-1]->Draw();
105     cD->cd(ip);
106     hNelDig[ip-1]->Draw();
107   }
108
109
110 }