Updated list of MUON libraries (Ivana)
[u/mrichter/AliRoot.git] / HMPID / Hqa.C
CommitLineData
00f49991 1void Hqa()
2{
3 gROOT->LoadMacro("Hdisp.C");
4 HitQA();
5}//Hqa()
6void HitQA()
7{
27311693 8 TFile *pFile = new TFile("Hqa.root","recreate");
9 AliHMPIDCluster::DoCorrSin(kFALSE);
10 AliHMPIDDigit::fSigmas=4;
11 Int_t nEvts=10000;
00f49991 12 TLegend *lQ=new TLegend(0.5,0.5,0.9,0.9);
13
14 TH1F *hQ7eV =new TH1F("hQ7eV" ,"" ,300,-50,2000); hQ7eV ->SetLineColor(kRed); lQ->AddEntry(hQ7eV ,"Ckov 7 eV"); hQ7eV->SetStats(0);
15 TH1F *hQ200eV=new TH1F("hQ200eV","" ,300,-50,2000); hQ200eV->SetLineColor(kBlack); lQ->AddEntry(hQ200eV,"mip 200 eV");
16 TH1F *hQ500eV=new TH1F("hQ500eV","" ,300,-50,2000); hQ500eV->SetLineColor(kCyan); lQ->AddEntry(hQ500eV,"mip 500 eV");
17 TH1F *hQ900eV=new TH1F("hQ900eV","" ,300,-50,2000); hQ900eV->SetLineColor(kGreen); lQ->AddEntry(hQ900eV,"mip 900 eV");
18
19 TH1F *hCluPerEvt=new TH1F("hCluPerEvt","# clusters per event",11,-0.5,10.5);
20 TH1F *hCluChi2 =new TH1F("hChi2","Chi2 ",1000,0,100);
21 TH1F *hCluFlg =new TH1F("hCluFlg","Cluster flag",14,-1.5,12.5); hCluFlg->SetFillColor(5);
22 TH1F *hCluRawSize= new TH1F("hCluRawSize","Raw cluster size ",100,0,100);
23
24 gStyle->SetOptStat(10);
25 TH2F *pCluMapSi1 =new TH2F("cluMapSi1","Size 1 map" ,1700,-10,160,1700,-10,160);
26 TH2F *pCluMapLo0 =new TH2F("cluMNoLo0","Loc Max 0 map" ,1700,-10,160,1700,-10,160);
27 TH2F *pCluMapLo1 =new TH2F("cluMapLo1","Loc Max 1 map" ,1700,-10,160,1700,-10,160);
28 TH2F *pCluMapUnf =new TH2F("cluMapUnf","Unfolded map" ,1700,-10,160,1700,-10,160);
29 TH2F *pCluMapEdg =new TH2F("cluMapEdg","On edge map" ,1700,-10,160,1700,-10,160);
30 TH2F *pCluMapCoG =new TH2F("cluMapCoG","CoG map" ,1700,-10,160,1700,-10,160);
31 TH2F *pCluMapEmp =new TH2F("cluMapEmp","undefined-empty" ,1700,-10,160,1700,-10,160);
32 TH2F *pCluMapNoLoc=new TH2F("cluMapNoLoc","no loc maxima" ,1700,-10,160,1700,-10,160);
33 TH2F *pCluMapAbn =new TH2F("cluMapAbn","abnormal fit" ,1700,-10,160,1700,-10,160);
34 TH2F *pCluMapNot =new TH2F("cluMapNot","Raw Clusters" ,1700,-10,160,1700,-10,160);
35 TH2F *pCluMapMax =new TH2F("cluMapMax","N. locs excceds" ,1700,-10,160,1700,-10,160);
36
27311693 37 TH1F *hHitCluDifX = new TH1F("hHitCluDifX" ,";entries;x_{Hit}-x_{Clu} [cm]" ,1000,-1,1); hHitCluDifX->Sumw2(); hHitCluDifX->SetFillColor(kYellow);
38// TH2F *hHitCluDifXv= new TH2F("hHitCluDifXv",";x_{Hit};x_{Hit}-x_{Clu} [cm]" ,500,-0.5,0.5,1000,-0.2,0.2);hHitCluDifXv->Sumw2();
39 TProfile *hHitCluDifXv= new TProfile("hHitCluDifXv",";x_{Hit};x_{Hit}-x_{Clu} [cm]" ,500,-0.5,0.5);
40 TH1F *hHitCluDifY = new TH1F("hHitCluDifY" ,";entries;y_{Hit}-y_{Clu} [cm]" ,1000,-1,1); hHitCluDifY->Sumw2(); hHitCluDifY->SetFillColor(kYellow);
41 TH2F *hHitCluDifXY= new TH2F("hHitCluDifXY",";x_{Hit}-x_{Clu};y_{Hit}-y_{Clu}",1000,-1,1,1000,-1,1);hHitCluDifXY->Sumw2();
00f49991 42 TH1F *hHitCluDifQ = new TH1F("hHitCluDifQ" ,";entries;(Q_{Clu}-Q_{Hit})/Q_{Hit}" ,200 ,-200,200); hHitCluDifQ->Sumw2(); hHitCluDifQ->SetFillColor(kYellow);
43
44 TH2F *hHitMap= new TH2F("hHitMap",";x_{Hit};y_{Hit}",1700,-10,160,1700,-10,160);
45
46 Float_t e200=200e-9,e500=500e-9,e900=900e-9,e7=7e-9;//predefined Eloss
47
48 AliHMPIDHit hit(0,0,kProton,0,0,0);
49 for(int i=0;i<5000;i++){
50 hQ200eV->Fill(hit.QdcTot(e200)); hQ500eV->Fill(hit.QdcTot(e500)); hQ900eV->Fill(hit.QdcTot(e900)); hQ7eV->Fill(hit.QdcTot(e7));
51 }
52 TClonesArray hits("AliHMPIDHit"); TClonesArray sdigs("AliHMPIDDigit");
53 TObjArray digs(7); for(Int_t i=0;i<7;i++) digs.AddAt(new TClonesArray("AliHMPIDDigit"),i);
54 TObjArray clus(7); for(Int_t i=0;i<7;i++) clus.AddAt(new TClonesArray("AliHMPIDCluster"),i);
55
56
57 for(Int_t iEvt=0;iEvt<nEvts;iEvt++){//events loop
27311693 58 if(iEvt%50==0)Printf("============> iEvt = %d ",iEvt);
00f49991 59
60 Int_t ch,pid; Float_t e,hitx,hity,hitq;
61// Int_t nHits=(type==999)?1:40;
62 Int_t nHits=1;
63 for(Int_t iHit=0;iHit<nHits;iHit++){//hits loop for the current event
64 switch(iHit){
27311693 65 case 0: ch=0;pid=kProton;e=e200;
66 hitx=AliHMPIDDigit::SizePadX()*(6+gRandom->Rndm());
67 hity=AliHMPIDDigit::SizePadY()*(6+gRandom->Rndm());break; //mip ramdomly distributed in one pad
00f49991 68 case 1: ch=0;pid=kProton;e=e200;hitx=0.4;hity=0.42;break; //mip in left-hand bottom coner of chamber 0
69 case 2: ch=0;pid=kProton;e=e200;hitx=0.4;hity=30 ;break; //mip on left edge of chamber 0
70 case 3: ch=0;pid=kProton;e=e200;hitx=40; hity=0.42;break; //mip on bottom edge of chamber 0
71 default: ch=gRandom->Rndm()*6; pid=(gRandom->Rndm()>0.9)? kProton:50000050;
72 if(pid==kProton)
73 e=gRandom->Rndm()*900e-9;
74 else
75 e=5.5e-9+3e-9*gRandom->Rndm();
76 hitx=gRandom->Rndm()*AliHMPIDDigit::SizeAllX(); hity=gRandom->Rndm()*AliHMPIDDigit::SizeAllY();break; //random hit
77 }
78 new(hits[iHit]) AliHMPIDHit(ch,e,pid,iHit,hitx,hity);
79 hitq=e;
80 }//hits loop
81
82 AliHMPIDv1::Hit2Sdi(&hits,&sdigs);
83 AliHMPIDDigitizer::DoNoise(kFALSE);
84 AliHMPIDDigitizer::Sdi2Dig(&sdigs,&digs);
27311693 85 AliHMPIDReconstructor::Dig2Clu(&digs,&clus,kTRUE);
00f49991 86
87// From here normal procedure for QA
88
89 for(Int_t iHit=0;iHit<hits.GetEntriesFast();iHit++) {
90 AliHMPIDHit *pHit = (AliHMPIDHit*)hits.UncheckedAt(iHit);
91 hHitMap->Fill(pHit->LorsX(),pHit->LorsY());
92 }
93 Int_t kMaxCh=(nHits==1)?0:AliHMPIDDigit::kMaxCh;
94 for(Int_t iCh=AliHMPIDDigit::kMinCh;iCh<=kMaxCh;iCh++){//chambers loop
95 TClonesArray *pDigs=(TClonesArray *)digs.UncheckedAt(iCh);
96 TClonesArray *pClus=(TClonesArray *)clus.UncheckedAt(iCh);
97
98// if(pClus->GetEntriesFast()>nHits) {hits.Print();}
99// if(pClus->GetEntriesFast()>nHits) {pDigs->Print();Printf("----S D I G I T S-------------");}
100// if(pClus->GetEntriesFast()>nHits) {sdigs.Print();Printf("-------------------------------");}
101 hCluPerEvt->Fill(pClus->GetEntriesFast());
102 for(Int_t iClu=0;iClu<pClus->GetEntriesFast();iClu++){//clusters loop
103 AliHMPIDCluster *pClu=(AliHMPIDCluster*)pClus->UncheckedAt(iClu);
104 Float_t clux=pClu->X(); Float_t cluy=pClu->Y(); Float_t cluq=pClu->Q();
105 hCluFlg->Fill(pClu->Status());
106 hCluChi2->Fill(pClu->Chi2());
107 hCluRawSize->Fill(pClu->Size());
108
109 switch(pClu->Status()){
110/*
111 case kFrm : status="formed " ;break;
112 case kUnf : status="unfolded (fit)" ;break;
113 case kCoG : status="coged " ;break;
114 case kLo1 : status="locmax 1 (fit)" ;break;
115 case kMax : status="exceeded (cog)" ;break;
116 case kNot : status="not done (cog)" ;break;
117 case kEmp : status="empty " ;break;
118 case kEdg : status="edge (fit)" ;break;
119 case kSi1 : status="size 1 (cog)" ;break;
120 case kNoLoc: status="no LocMax(fit)" ;break;
121 case kAbn : status="Abnormal fit " ;break;
122*/
27311693 123 case AliHMPIDCluster::kSi1: pCluMapSi1->Fill(clux,cluy); break;
124 case AliHMPIDCluster::kLo1: pCluMapLo1->Fill(clux,cluy); break;
125 case AliHMPIDCluster::kUnf: pCluMapUnf->Fill(clux,cluy); break;
126 case AliHMPIDCluster::kMax: pCluMapMax->Fill(clux,cluy); break;
127 case AliHMPIDCluster::kEdg: pCluMapEdg->Fill(clux,cluy); break;
128 case AliHMPIDCluster::kCoG: pCluMapCoG->Fill(clux,cluy); break;
129 case AliHMPIDCluster::kNoLoc: pCluMapNoLoc->Fill(clux,cluy);break;
130 case AliHMPIDCluster::kAbn: pCluMapAbn->Fill(clux,cluy); break;
131 case AliHMPIDCluster::kNot: pCluMapNot->Fill(clux,cluy); break;
00f49991 132 default: pCluMapEmp->Fill(clux,cluy); break;
133 }
134
135 hHitCluDifX->Fill(hitx-clux); hHitCluDifY->Fill(hity-cluy); hHitCluDifXY->Fill(hitx-clux,hity-cluy); hHitCluDifQ->Fill(100*(cluq-hitq)/hitq);
27311693 136 // distorsion due to feedback photons
137 Int_t pc,px,py;
138 AliHMPIDDigit::Lors2Pad(hitx,hity,pc,px,py);
139 Float_t padCenterX = AliHMPIDDigit::LorsX(pc,px);
140 if(pClu->Size()>1)hHitCluDifXv->Fill(hitx-padCenterX,(hitx-clux));
141 //
00f49991 142 }//clusters loop
143 }//chambers loop
144
145 hits.Delete(); sdigs.Delete(); for(int i = 0;i<7;i++){((TClonesArray*)digs.At(i))->Delete();((TClonesArray*)clus.At(i))->Delete();}
146 }//events loop
147
148 gStyle->SetPalette(1);
149 TCanvas *pC2=new TCanvas("Digit canvas","Digit canvas",1280,800); pC2->Divide(3,3);
150 pC2->cd(1);hHitCluDifX->Draw("hist");
151 pC2->cd(2);hHitCluDifY->Draw("hist");
27311693 152// pC2->cd(3);hHitCluDifXY->Draw("colz");
153 // Draw CorrSin
154 AliHMPIDCluster c;
155 TPolyLine *pLine = new TPolyLine(500);
156 for(Int_t i=0;i<500;i++) {
157 Double_t x = 0 + i*AliHMPIDDigit::SizePadX()/500.;
158 c.SetX(x);c.SetY(0);c.CorrSin();
159 pLine->SetPoint(i,x-0.5*AliHMPIDDigit::SizePadX(),c.X()-x);
160 }
161 pC2->cd(3);hHitCluDifXv->Draw("colz");pLine->Draw("C");
162 //
00f49991 163 pC2->cd(4);hHitCluDifQ->Draw("hist");
164 pC2->cd(5);gPad->SetLogy(1);hCluFlg->Draw();
165 pC2->cd(6);hCluChi2->Draw();
166 pC2->cd(7); hCluRawSize->Draw();
167 pC2->cd(8); hCluPerEvt->Draw();
168 pC2->cd(9); hQ7eV->Draw(); hQ200eV->Draw("same"); hQ500eV->Draw("same"); hQ900eV->Draw("same"); lQ->Draw();
169 TCanvas *pC1=new TCanvas("ClusterMaps","Cluster maps",1280,800); pC1->Divide(3,3);
170 pC1->cd(1); pCluMapSi1->Draw(); DrawCh(-1);
171 pC1->cd(2); pCluMapLo1->Draw(); DrawCh(-1);
172 pC1->cd(3); pCluMapUnf->Draw(); DrawCh(-1);
173 pC1->cd(4); hHitMap->Draw(); DrawCh(-1);
174 pC1->cd(5); pCluMapMax->Draw(); DrawCh(-1);
175 pC1->cd(6); pCluMapEdg->Draw(); DrawCh(-1);
176 pC1->cd(7); pCluMapNot->Draw(); DrawCh(-1);
177 pC1->cd(8); pCluMapNoLoc->Draw();DrawCh(-1);
178 pC1->cd(9); pCluMapEmp->Draw(); DrawCh(-1);
179
180 pC1->SaveAs("$HOME/HitMaps.png"); //?????
181 pC2->SaveAs("$HOME/HitCluDif.gif");
27311693 182
183 pFile->Write();
00f49991 184 Printf("Digits - raw -digits conversion...");
185
186 AliHMPIDDigit d1,d2; Int_t ddl,r,d,a;UInt_t w32;
187
188 for(Int_t iCh=AliHMPIDDigit::kMinCh;iCh<=AliHMPIDDigit::kMaxCh;iCh++)
189 for(Int_t iPc=AliHMPIDDigit::kMinPc;iPc<=AliHMPIDDigit::kMaxPc;iPc++)
190 for(Int_t iPx=AliHMPIDDigit::kMinPx;iPx<=AliHMPIDDigit::kMaxPx;iPx++)
191 for(Int_t iPy=AliHMPIDDigit::kMinPy;iPy<=AliHMPIDDigit::kMaxPy;iPy++){
192 d1.Set(iCh,iPc,iPx,iPy,3040); //set digit
193 d1.Raw(w32,ddl,r,d,a); //get raw word for this digit
194 d2.Raw(w32,ddl); //set another digit from that raw word
195 if(d1.Compare(&d2)) {d1.Print(); d2.Print(); Printf("");}//compare them
196 }
197 Printf("OK");
198}//tst()