More checks
[u/mrichter/AliRoot.git] / HMPID / Hqa.C
1 #if !defined(__CINT__) || defined(__MAKECINT__)
2 #include <TSystem.h>
3 #include <TFile.h>
4 #include <TStyle.h>
5 #include <TTree.h>
6 #include <TClonesArray.h>
7 #include <TObjArray.h>
8 #include <TF1.h>
9 #include <TH1F.h>
10 #include <TH2F.h>
11 #include <TProfile.h>
12 #include <TCanvas.h>
13
14 #include <AliHMPIDParam.h>
15 #include <AliHMPIDHit.h>
16 #include <AliHMPIDCluster.h>
17 #include <AliHMPIDDigit.h>
18 #endif
19
20 Int_t nEntries = 0;
21
22 TObjArray *CreateContainer(const char *classname,TTree *pTree);
23 void Hits(Int_t mode,TTree *pTree=0x0);
24 void Digs(Int_t mode, TTree *pTree=0x0);
25 void Clus(Int_t mode, TTree *pTree=0x0);
26 void Summary();
27   
28 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
29 TObjArray *CreateContainer(const char *classname,TTree *pTree)
30 {
31   TObjArray *pOA=new TObjArray(AliHMPIDParam::kMaxCh+1);
32   for(Int_t iCh=AliHMPIDParam::kMinCh;iCh<=AliHMPIDParam::kMaxCh;iCh++){
33     TClonesArray *pCA=new TClonesArray(classname);
34     pOA->AddAt(pCA,iCh);    
35   }
36   
37   pOA->SetOwner(kTRUE);  
38   
39   for(Int_t iCh=AliHMPIDParam::kMinCh;iCh<=AliHMPIDParam::kMaxCh;iCh++){
40     pTree->SetBranchAddress(Form("HMPID%i",iCh),&(*pOA)[iCh]);
41   }
42   return pOA;
43 }
44 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
45 TH1F *hHitQdc,*hHitQdcCh[AliHMPIDParam::kMaxCh+1]; TH2F *hHitMap[AliHMPIDParam::kMaxCh+1];  
46 Double_t fHitMean[AliHMPIDParam::kMaxCh+1],fHitErr[AliHMPIDParam::kMaxCh+1];
47 void Hits(Int_t mode,TTree *pTree)
48 {
49   switch(mode){
50     case 1:
51       hHitQdc=new TH1F("HitQdc","Hit Qdc all chamber;QDC",4000,0,4000);
52   for(Int_t iCh=AliHMPIDParam::kMinCh;iCh<=AliHMPIDParam::kMaxCh;iCh++){
53         hHitMap[iCh]=new TH2F(Form("HitMap%i",iCh),Form("Ch%i;x_{Hit};y_{Hit}",iCh),160,0,160,160,0,160);      
54         hHitQdcCh[iCh]=new TH1F(Form("HMPID%i",iCh),Form("Charge for HMPID%i",iCh),4000,0,4000);
55       }
56       break;
57     case 2:
58       if(pTree==0) return;
59       TClonesArray *pHits=new TClonesArray("AliHMPIDHit");  pTree->SetBranchAddress("HMPID",&pHits);  
60       for(Int_t iEnt=0;iEnt<pTree->GetEntriesFast();iEnt++){//entries loop
61         pTree->GetEntry(iEnt);
62         for(Int_t iHit=0;iHit<pHits->GetEntriesFast();iHit++){//hits loop
63           AliHMPIDHit *pHit = (AliHMPIDHit*)pHits->UncheckedAt(iHit);
64           hHitMap[pHit->Ch()]->Fill(pHit->LorsX(),pHit->LorsY());
65           hHitQdc->Fill(pHit->Q());
66           hHitQdcCh[pHit->Ch()]->Fill(pHit->Q());
67         }//hits loop      
68       }//entries loop
69       delete pHits;
70       break;
71     case 3:
72       TCanvas *c1=new TCanvas("HitCan","Hits",1280,800); c1->Divide(3,3);
73   
74       for(Int_t iCh=AliHMPIDParam::kMinCh;iCh<=AliHMPIDParam::kMaxCh;iCh++){
75         if(iCh==6) c1->cd(1); if(iCh==5) c1->cd(2);
76         if(iCh==4) c1->cd(4); if(iCh==3) c1->cd(5); if(iCh==2) c1->cd(6);
77                               if(iCh==1) c1->cd(8); if(iCh==0) c1->cd(9);
78         gStyle->SetPalette(1);      
79         hHitMap[iCh]->Draw("colz");
80       }  
81       for(Int_t iCh=AliHMPIDParam::kMinCh;iCh<=AliHMPIDParam::kMaxCh;iCh++){
82         fHitMean[iCh] = 0;
83         fHitErr[iCh]  = 0;
84         if((Int_t)hHitQdcCh[iCh]->GetEntries()<nEntries) continue;
85         c1->cd(3);hHitQdcCh[iCh]->Fit("expo","Q");
86         TF1 *funcfit = (TF1*)hHitQdcCh[iCh]->FindObject("expo");
87         fHitMean[iCh] = funcfit->GetParameter(1);
88         fHitErr[iCh]  = funcfit->GetParError(1);
89       }
90       TPad *pad = (TPad*)c1->cd(3); hHitQdc->SetFillColor(5); pad->SetLogy(); hHitQdc->Fit("expo","Q");
91       break;
92   }
93 }//Hits()
94 //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
95 TH1F     *hDigQ;
96 TProfile *hDigHighQ,*hDigChEvt;
97 void Digs(Int_t mode, TTree *pTree)
98 {
99   switch(mode){
100     case 1:
101       hDigHighQ=new TProfile("hDigHighQ","Highest charge in chamber  ",AliHMPIDParam::kMaxCh+1,AliHMPIDParam::kMinCh,AliHMPIDParam::kMaxCh+1);
102       hDigChEvt=new TProfile("hDigChEvt","Chamber occupancy per event",AliHMPIDParam::kMaxCh+1,AliHMPIDParam::kMinCh,AliHMPIDParam::kMaxCh+1);
103       hDigQ    =new TH1F("hDigQ        ","Charge of digits (ADC)     ",3000,0,3000);
104       break;
105     case 2:
106       if(pTree==0) return;
107       TObjArray *pLst=CreateContainer("AliHMPIDDigit",pTree); pTree->GetEntry(0);
108       for(Int_t iCh=AliHMPIDParam::kMinCh;iCh<=AliHMPIDParam::kMaxCh;iCh++){//chambers loop
109         TClonesArray *pDigs=(TClonesArray *)pLst->UncheckedAt(iCh);
110         hDigChEvt->Fill(iCh,pDigs->GetEntriesFast()/(48.*80.*6.)*100.);
111         Double_t highQ = 0;
112         for(Int_t iDig=0;iDig<pDigs->GetEntriesFast();iDig++){//digits loop
113           AliHMPIDDigit *pDig=(AliHMPIDDigit*)pDigs->UncheckedAt(iDig);
114           hDigQ->Fill(pDig->Q());
115           if(pDig->Q()>highQ) highQ = pDig->Q();
116         }
117         hDigHighQ->Fill(iCh,highQ);
118       }
119       delete pLst;
120       break;
121     case 3:
122       TCanvas *c1=new TCanvas("DigQa","Digit Check",1280,800); c1->Divide(2,2);
123       gStyle->SetOptFit(1);
124       TPad *pad = (TPad*)c1->cd(1); pad->SetLogy(); hDigQ->Fit("expo","QN"); 
125       c1->cd(2); hDigHighQ->Draw();
126       c1->cd(3); hDigChEvt->Draw();
127       break;
128   }//switch
129 }//Dig()
130 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
131 TH1F *hCluEvt,*hCluChi2,*hCluFlg,*hCluSize,*hCluQ;
132 void Clus(Int_t mode, TTree *pTree)
133 {
134   switch(mode){
135     case 1:
136       hCluEvt=new TH1F("CluPerEvt","Cluster multiplicity"   ,100,0,100);
137       hCluChi2  =new TH1F("CluChi2"  ,"Chi2 "               ,1000,0,100);
138       hCluFlg   =new TH1F("CluFlg"   ,"Cluster flag"        ,14,-1.5,12.5);                       hCluFlg->SetFillColor(5);
139       hCluSize  =new TH1F("CluSize"  ,"Cluster size        ",100,0,100);
140       hCluQ     =new TH1F("CluQ"     ,"Cluster charge (ADC)",1000,0,5000);
141       break;
142     case 2:      
143       if(pTree==0) return;
144
145       TObjArray *pLst=CreateContainer("AliHMPIDCluster",pTree); pTree->GetEntry(0);
146       for(Int_t iCh=AliHMPIDParam::kMinCh;iCh<=AliHMPIDParam::kMaxCh;iCh++){//chambers loop
147         TClonesArray *pClus=(TClonesArray *)pLst->UncheckedAt(iCh);
148         hCluEvt->Fill(pClus->GetEntriesFast());
149         for(Int_t iClu=0;iClu<pClus->GetEntriesFast();iClu++){//clusters loop
150           AliHMPIDCluster *pClu=(AliHMPIDCluster*)pClus->UncheckedAt(iClu);
151           hCluFlg->Fill(pClu->Status());
152           hCluChi2->Fill(pClu->Chi2());
153           hCluSize->Fill(pClu->Size());
154           hCluQ->Fill(pClu->Q());
155         }
156       }
157       delete pLst;           
158       break;
159     case 3:
160       TCanvas *c1=new TCanvas("CluComCan","Clusters in common",1280,800); c1->Divide(3,3);
161       c1->cd(1); hCluEvt->SetFillColor(5);      hCluEvt->Draw();
162       c1->cd(2); hCluChi2->SetFillColor(5);     hCluChi2->Draw(); 
163       c1->cd(3); hCluFlg->SetFillColor(5);      hCluFlg->Draw(); 
164       c1->cd(4); hCluSize->SetFillColor(5);     hCluSize->Draw(); 
165       TPad *pad = (TPad*)c1->cd(5); hCluQ->SetFillColor(5); pad->SetLogy(); hCluQ->Draw(); 
166       break;
167   }//switch
168 }//Clus()
169 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
170 void Hqa()
171 {
172   TFile *fh=0; if(gSystem->IsFileInIncludePath("HMPID.Hits.root"))      fh=TFile::Open("HMPID.Hits.root"     ,"read");if(fh->IsZombie()) fh=0;
173   TFile *fd=0; if(gSystem->IsFileInIncludePath("HMPID.Digits.root"))    fd=TFile::Open("HMPID.Digits.root"   ,"read");if(fd->IsZombie()) fd=0;
174   TFile *fc=0; if(gSystem->IsFileInIncludePath("HMPID.RecPoints.root")) fc=TFile::Open("HMPID.RecPoints.root","read");if(fc->IsZombie()) fc=0;
175   if(fh==0 && fd==0 && fc==0){Printf("Nothing to do!"); return;}
176   if(fh) Hits(1); if(fc) Clus(1);  if(fd) Digs(1);//book
177   Int_t iEvt=0;
178   while(1){
179     TTree *th=0; if(fh) th=(TTree*)fh->Get(Form("Event%i/TreeH",iEvt));
180     TTree *td=0; if(fd) td=(TTree*)fd->Get(Form("Event%i/TreeD",iEvt));
181     TTree *tc=0; if(fc) tc=(TTree*)fc->Get(Form("Event%i/TreeR",iEvt));
182     
183     Hits(2,th); Digs(2,td); Clus(2,tc); //fill
184     if(th==0 && td==0 && tc==0) break;
185     iEvt++;
186     if(!(iEvt%50)) Printf("Event %i processed",iEvt);
187   }
188   
189   nEntries = iEvt;
190   
191   if(fd) Clus(3);//plot everything
192   if(fc) Digs(3); 
193   if(fh) Hits(3);
194   Summary();
195 }
196 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
197 void Summary()
198 {
199   //info for hits...
200   for(Int_t iCh=AliHMPIDParam::kMinCh;iCh<=AliHMPIDParam::kMaxCh;iCh++){
201     Printf(" #################### Summary for HMPID %i#################### ",iCh);
202     //info for hits...
203     Printf("-----Summary of Hits----- ");
204     if(fHitMean[iCh]==0) {
205       Printf("gain %5.2f +/- %5.2f",fHitMean[iCh],fHitErr[iCh]);
206     } else {   
207       Double_t gain = 1./TMath::Abs(fHitMean[iCh]);
208       Double_t errgain = gain*gain*fHitErr[iCh];
209       Printf("gain %5.2f +/- %5.2f",gain,errgain);
210     }
211     //info for digits...
212     Printf("-----Summary of Digits-----");
213     Printf(" Chamber %d with occupancy (%) %5.2f +/- %5.2f",iCh,hDigChEvt->GetBinContent(iCh+1),hDigChEvt->GetBinError(iCh+1));
214     Printf(" Chamber %d with higest Q (ADC) %7.0f +/- %7.0f",iCh,hDigHighQ->GetBinContent(iCh+1),hDigHighQ->GetBinError(iCh+1));
215   }
216 }