new histos
[u/mrichter/AliRoot.git] / HMPID / Hqa.C
1 #include <TSystem.h>
2 #include <TFile.h>
3 #include <TTree.h>
4 #include <TClonesArray.h>
5 #include <TObjArray.h>
6 #include <TH2F.h>
7 #include <TCanvas.h>
8
9 #include <AliHMPIDHit.h>
10 #include <AliHMPIDCluster.h>
11
12
13
14 TH1F *hHitQdc; TH2F *hHitMap[7];  
15
16
17
18 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
19
20 void BookHits()
21 {
22   hHitQdc=new TH1F("HitQdc","Hit Qdc all chamber;QDC",500,0,4000);
23   for(Int_t iCh=0;iCh<7;iCh++)
24     hHitMap[iCh]=new TH2F(Form("HitMap%i",iCh),Form("Ch%i;x_{Hit};y_{Hit}",iCh),1700,-10,160,1700,-10,160);  
25 }
26 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
27 TObjArray *CreateContainer(const char *classname,TTree *pTree)
28 {
29   TObjArray *pOA=new TObjArray(7); pOA->SetOwner();
30   for(Int_t iCh=AliHMPIDDigit::kMinCh;iCh<=AliHMPIDDigit::kMaxCh;iCh++){
31     TClonesArray *pCA=new TClonesArray(classname);
32     pOA->AddAt(pCA,iCh);    
33     pTree->SetBranchAddress(Form("HMPID%i",iCh),&pCA);
34   }
35   pTree->GetEntry(0);
36   return pOA;
37 }
38 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
39
40
41 TH1F *hCluEvt,*hCluChi2,*hCluFlg,*hCluSize;
42
43 void Clus(Int_t mode, TTree *pTree=0x0)
44 {
45   switch(mode){
46     case 1:
47       hCluEvt=new TH1F("CluPerEvt","# clusters per event",11,-0.5,10.5);
48       hCluChi2  =new TH1F("CluChi2"  ,"Chi2 "               ,1000,0,100);
49       hCluFlg   =new TH1F("CluFlg"   ,"Cluster flag"        ,14,-1.5,12.5);                       hCluFlg->SetFillColor(5);
50       hCluSize  =new TH1F("CluSize"  ,"Raw cluster size    ",100,0,100);
51       return;
52     case 2:      
53       if(pTree==0) return;
54       
55       TObjArray *pLst=CreateContainer("AliHMPIDCluster",pTree); 
56       for(Int_t iCh=AliHMPIDDigit::kMinCh;iCh<=AliHMPIDDigit::kMaxCh;iCh++){//chambers loop
57         TClonesArray *pClus=(TClonesArray *)pLst->UncheckedAt(iCh);
58       
59         hCluEvt->Fill(pClus->GetEntriesFast());
60         for(Int_t iClu=0;iClu<pClus->GetEntriesFast();iClu++){//clusters loop
61           AliHMPIDCluster *pClu=(AliHMPIDCluster*)pClus->UncheckedAt(iClu);
62           hCluFlg->Fill(pClu->Status());
63           hCluChi2->Fill(pClu->Chi2());
64           hCluSize->Fill(pClu->Size());
65         }
66       }
67       delete pLst;           
68       return;  
69     case 3:
70       TCanvas *c1=new TCanvas("CluComCan","Clusters in common",1280,800); c1->Divide(3,3);
71       c1->cd(1); hCluEvt->Draw();  
72       c1->cd(2); hCluChi2->Draw(); 
73       c1->cd(3); hCluFlg->Draw(); 
74       c1->cd(4); hCluSize->Draw(); 
75       return;
76   }//switch
77 }//Clus()
78 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
79
80 void FillHits(TTree *pHitTree)
81 {
82   if(pHitTree==0) return;
83   TClonesArray *pHitLst=new TClonesArray("AliHMPIDHit");
84   pHitTree->SetBranchAddress("HMPID",&pHitLst);
85   
86   for(Int_t iEnt=0;iEnt<pHitTree->GetEntriesFast();iEnt++){
87     pHitTree->GetEntry(iEnt);
88     for(Int_t iHit=0;iHit<pHitLst->GetEntriesFast();iHit++){
89       AliHMPIDHit *pHit = (AliHMPIDHit*)pHitLst->UncheckedAt(iHit);
90       hHitMap[pHit->Ch()]->Fill(pHit->LorsX(),pHit->LorsY());
91       hHitQdc->Fill(pHit->Q());
92     }      
93   }
94   delete pHitLst;
95 }    
96 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
97
98 void PlotHits()
99 {
100   TCanvas *pHitCan=new TCanvas("HitCan","Hits",1280,800); pHitCan->Divide(3,3);
101   
102   for(Int_t iCh=0;iCh<7;iCh++){
103     if(iCh==6) pHitCan->cd(1); if(iCh==5) pHitCan->cd(2);
104     if(iCh==4) pHitCan->cd(4); if(iCh==3) pHitCan->cd(5); if(iCh==2) pHitCan->cd(6);
105                                if(iCh==1) pHitCan->cd(8); if(iCh==0) pHitCan->cd(9);
106     hHitMap[iCh]->Draw();
107   }  
108   pHitCan->cd(3); gPad->SetLogy(); hHitQdc->Draw();
109 }
110
111 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
112
113 void Hqa()
114 {
115
116   TFile *fh=0; if(gSystem->IsFileInIncludePath("HMPID.Hits.root"))      fh=TFile::Open("HMPID.Hits.root"     ,"read");if(fh->IsZombie()) fh=0;
117   TFile *fd=0; if(gSystem->IsFileInIncludePath("HMPID.Digits.root"))    fd=TFile::Open("HMPID.Digits.root"   ,"read");if(fd->IsZombie()) fd=0;
118   TFile *fc=0; if(gSystem->IsFileInIncludePath("HMPID.RecPoints.root")) fc=TFile::Open("HMPID.RecPoints.root","read");if(fc->IsZombie()) fc=0;
119   
120   
121   if(fh==0 && fd==0 && fc==0){Printf("Nothing to do!"); return;}
122   
123   
124   if(fh) BookHits();
125   if(fc) Clus(1); //book
126   
127   Int_t iEvt=0;
128   while(1){
129     TTree *th=0; if(fh) th=(TTree*)fh->Get(Form("Event%i/TreeH",iEvt));
130     TTree *td=0; if(fd) td=(TTree*)fd->Get(Form("Event%i/TreeD",iEvt));
131     TTree *tc=0; if(fc) tc=(TTree*)fc->Get(Form("Event%i/TreeR",iEvt));
132     
133     FillHits(th);
134     Clus(2,tc);
135     if(th==0 && td==0 && tc==0) break;
136     iEvt++;
137     Printf("Event %i processed",iEvt);
138   }
139   
140   if(fh) PlotHits();
141   if(fc) Clus(3); //plot
142 }
143
144 /*
145
146 void BlahQA()
147 {
148   TFile *pFile = new TFile("Hqa.root","recreate");
149   AliHMPIDCluster::DoCorrSin(kFALSE);
150   AliHMPIDDigit::fSigmas=4;
151   Int_t nEvts=10000;
152   TLegend *lQ=new TLegend(0.5,0.5,0.9,0.9);
153
154   TH1F *hQ7eV  =new TH1F("hQ7eV"  ,"" ,300,-50,2000);     hQ7eV  ->SetLineColor(kRed);      lQ->AddEntry(hQ7eV  ,"Ckov 7 eV");  hQ7eV->SetStats(0);
155   TH1F *hQ200eV=new TH1F("hQ200eV","" ,300,-50,2000);     hQ200eV->SetLineColor(kBlack);    lQ->AddEntry(hQ200eV,"mip 200 eV");
156   TH1F *hQ500eV=new TH1F("hQ500eV","" ,300,-50,2000);     hQ500eV->SetLineColor(kCyan);     lQ->AddEntry(hQ500eV,"mip 500 eV");
157   TH1F *hQ900eV=new TH1F("hQ900eV","" ,300,-50,2000);     hQ900eV->SetLineColor(kGreen);    lQ->AddEntry(hQ900eV,"mip 900 eV");
158   
159   
160   gStyle->SetOptStat(10);
161   TH2F *pCluMapSi1  =new TH2F("cluMapSi1","Size 1 map"       ,1700,-10,160,1700,-10,160);
162   TH2F *pCluMapLo0  =new TH2F("cluMNoLo0","Loc Max 0 map"    ,1700,-10,160,1700,-10,160);     
163   TH2F *pCluMapLo1  =new TH2F("cluMapLo1","Loc Max 1 map"    ,1700,-10,160,1700,-10,160);        
164   TH2F *pCluMapUnf  =new TH2F("cluMapUnf","Unfolded map"     ,1700,-10,160,1700,-10,160);         
165   TH2F *pCluMapEdg  =new TH2F("cluMapEdg","On edge map"      ,1700,-10,160,1700,-10,160);         
166   TH2F *pCluMapCoG  =new TH2F("cluMapCoG","CoG  map"         ,1700,-10,160,1700,-10,160);            
167   TH2F *pCluMapEmp  =new TH2F("cluMapEmp","undefined-empty"  ,1700,-10,160,1700,-10,160);      
168   TH2F *pCluMapNoLoc=new TH2F("cluMapNoLoc","no loc maxima"  ,1700,-10,160,1700,-10,160);      
169   TH2F *pCluMapAbn  =new TH2F("cluMapAbn","abnormal fit"     ,1700,-10,160,1700,-10,160);      
170   TH2F *pCluMapNot  =new TH2F("cluMapNot","Raw Clusters"     ,1700,-10,160,1700,-10,160);      
171   TH2F *pCluMapMax  =new TH2F("cluMapMax","N. locs excceds"  ,1700,-10,160,1700,-10,160);      
172
173   TH1F *hHitCluDifX = new TH1F("hHitCluDifX" ,";entries;x_{Hit}-x_{Clu} [cm]"   ,1000,-1,1);          hHitCluDifX->Sumw2();    hHitCluDifX->SetFillColor(kYellow);
174 //  TH2F *hHitCluDifXv= new TH2F("hHitCluDifXv",";x_{Hit};x_{Hit}-x_{Clu} [cm]"   ,500,-0.5,0.5,1000,-0.2,0.2);hHitCluDifXv->Sumw2();
175   TProfile *hHitCluDifXv= new TProfile("hHitCluDifXv",";x_{Hit};x_{Hit}-x_{Clu} [cm]"   ,500,-0.5,0.5);
176   TH1F *hHitCluDifY = new TH1F("hHitCluDifY" ,";entries;y_{Hit}-y_{Clu} [cm]"   ,1000,-1,1);          hHitCluDifY->Sumw2();    hHitCluDifY->SetFillColor(kYellow);
177   TH2F *hHitCluDifXY= new TH2F("hHitCluDifXY",";x_{Hit}-x_{Clu};y_{Hit}-y_{Clu}",1000,-1,1,1000,-1,1);hHitCluDifXY->Sumw2();  
178   TH1F *hHitCluDifQ = new TH1F("hHitCluDifQ" ,";entries;(Q_{Clu}-Q_{Hit})/Q_{Hit}" ,200 ,-200,200);   hHitCluDifQ->Sumw2();    hHitCluDifQ->SetFillColor(kYellow);
179   
180  
181   Float_t e200=200e-9,e500=500e-9,e900=900e-9,e7=7e-9;//predefined  Eloss
182   
183   AliHMPIDHit hit(0,0,kProton,0,0,0);
184   for(int i=0;i<5000;i++){
185     hQ200eV->Fill(hit.QdcTot(e200));  hQ500eV->Fill(hit.QdcTot(e500));   hQ900eV->Fill(hit.QdcTot(e900));  hQ7eV->Fill(hit.QdcTot(e7));
186   }  
187   TClonesArray hits("AliHMPIDHit");  TClonesArray sdigs("AliHMPIDDigit");
188   TObjArray digs(7); for(Int_t i=0;i<7;i++) digs.AddAt(new TClonesArray("AliHMPIDDigit"),i);
189   TObjArray clus(7); for(Int_t i=0;i<7;i++) clus.AddAt(new TClonesArray("AliHMPIDCluster"),i);
190   
191     
192   for(Int_t iEvt=0;iEvt<nEvts;iEvt++){//events loop
193     if(iEvt%50==0)Printf("============> iEvt = %d ",iEvt);
194     
195     Int_t ch,pid; Float_t e,hitx,hity,hitq;
196 //    Int_t nHits=(type==999)?1:40;
197     Int_t nHits=1;
198     for(Int_t iHit=0;iHit<nHits;iHit++){//hits loop for the current event
199       switch(iHit){
200         case 0:  ch=0;pid=kProton;e=e200;
201                  hitx=AliHMPIDDigit::SizePadX()*(6+gRandom->Rndm());
202                  hity=AliHMPIDDigit::SizePadY()*(6+gRandom->Rndm());break; //mip ramdomly distributed in one pad
203         case 1:  ch=0;pid=kProton;e=e200;hitx=0.4;hity=0.42;break; //mip in left-hand bottom coner of chamber 0
204         case 2:  ch=0;pid=kProton;e=e200;hitx=0.4;hity=30  ;break; //mip on left edge of chamber 0
205         case 3:  ch=0;pid=kProton;e=e200;hitx=40; hity=0.42;break; //mip on bottom edge of chamber 0
206         default: ch=gRandom->Rndm()*6; pid=(gRandom->Rndm()>0.9)? kProton:50000050;
207                   if(pid==kProton) 
208                     e=gRandom->Rndm()*900e-9; 
209                   else
210                     e=5.5e-9+3e-9*gRandom->Rndm();
211                   hitx=gRandom->Rndm()*AliHMPIDDigit::SizeAllX(); hity=gRandom->Rndm()*AliHMPIDDigit::SizeAllY();break; //random hit
212       }
213       new(hits[iHit]) AliHMPIDHit(ch,e,pid,iHit,hitx,hity);                          
214       hitq=e;
215     }//hits loop
216     
217     AliHMPIDv1::Hit2Sdi(&hits,&sdigs);
218     AliHMPIDDigitizer::DoNoise(kFALSE);
219     AliHMPIDDigitizer::Sdi2Dig(&sdigs,&digs);
220     AliHMPIDReconstructor::Dig2Clu(&digs,&clus,kTRUE);
221
222 // From here normal procedure for QA
223
224     Int_t kMaxCh=(nHits==1)?0:AliHMPIDDigit::kMaxCh;
225     for(Int_t iCh=AliHMPIDDigit::kMinCh;iCh<=kMaxCh;iCh++){//chambers loop
226       TClonesArray *pDigs=(TClonesArray *)digs.UncheckedAt(iCh);
227       TClonesArray *pClus=(TClonesArray *)clus.UncheckedAt(iCh);
228       
229 //      if(pClus->GetEntriesFast()>nHits) {hits.Print();}
230 //      if(pClus->GetEntriesFast()>nHits) {pDigs->Print();Printf("----S D I G I T S-------------");}
231 //      if(pClus->GetEntriesFast()>nHits) {sdigs.Print();Printf("-------------------------------");}
232       hCluPerEvt->Fill(pClus->GetEntriesFast());
233       for(Int_t iClu=0;iClu<pClus->GetEntriesFast();iClu++){//clusters loop
234         AliHMPIDCluster *pClu=(AliHMPIDCluster*)pClus->UncheckedAt(iClu);
235         Float_t clux=pClu->X(); Float_t cluy=pClu->Y(); Float_t cluq=pClu->Q();
236         hCluFlg->Fill(pClu->Status());
237         hCluChi2->Fill(pClu->Chi2());
238         hCluRawSize->Fill(pClu->Size());
239          
240         switch(pClu->Status()){
241             case AliHMPIDCluster::kSi1:   pCluMapSi1->Fill(clux,cluy);  break;
242             case AliHMPIDCluster::kLo1:   pCluMapLo1->Fill(clux,cluy);  break;
243             case AliHMPIDCluster::kUnf:   pCluMapUnf->Fill(clux,cluy);  break; 
244             case AliHMPIDCluster::kMax:   pCluMapMax->Fill(clux,cluy);  break;
245             case AliHMPIDCluster::kEdg:   pCluMapEdg->Fill(clux,cluy);  break;
246             case AliHMPIDCluster::kCoG:   pCluMapCoG->Fill(clux,cluy);  break;
247             case AliHMPIDCluster::kNoLoc: pCluMapNoLoc->Fill(clux,cluy);break;
248             case AliHMPIDCluster::kAbn:   pCluMapAbn->Fill(clux,cluy);  break;
249             case AliHMPIDCluster::kNot:   pCluMapNot->Fill(clux,cluy);  break;
250             default:     pCluMapEmp->Fill(clux,cluy); break;            
251         }
252         
253         hHitCluDifX->Fill(hitx-clux); hHitCluDifY->Fill(hity-cluy); hHitCluDifXY->Fill(hitx-clux,hity-cluy); hHitCluDifQ->Fill(100*(cluq-hitq)/hitq);
254         // distorsion due to feedback photons
255         Int_t pc,px,py;
256         AliHMPIDDigit::Lors2Pad(hitx,hity,pc,px,py);
257         Float_t padCenterX = AliHMPIDDigit::LorsX(pc,px);
258         if(pClu->Size()>1)hHitCluDifXv->Fill(hitx-padCenterX,(hitx-clux));        
259         //
260       }//clusters loop
261     }//chambers loop
262       
263     hits.Delete();  sdigs.Delete();  for(int i = 0;i<7;i++){((TClonesArray*)digs.At(i))->Delete();((TClonesArray*)clus.At(i))->Delete();}
264   }//events loop      
265 }
266
267
268
269 */