- Fix problem with number of events to treat;
[u/mrichter/AliRoot.git] / HMPID / Hdisp.C
CommitLineData
8282f9d0 1TCanvas *pAll=0;
2AliRunLoader *gAL=0; AliLoader *gHL=0; AliESD *gEsd=0; TTree *gEsdTr=0; AliHMPID *gH=0;
3Int_t gEvt=0; Int_t gMaxEvt=0;
bd73b83f 4TObjArray *pNmean;
8282f9d0 5//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6void Hdisp()
7{//display events from files if any in current directory or simulated events
50e69536 8 pAll=new TCanvas("all","",1300,900); pAll->Divide(3,3,0,0);
9// pAll->ToggleEditor();
8282f9d0 10 pAll->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)",0,"","DoZoom(Int_t,Int_t,Int_t,TObject*)");
bd73b83f 11
12 OpenCalib();
8282f9d0 13 if(gSystem->IsFileInIncludePath("galice.root")){// tries to open session
14 if(gAlice) delete gAlice; //in case we execute this in aliroot delete default AliRun object
15 gAL=AliRunLoader::Open(); //try to open galice.root from current dir
16 gAL->LoadgAlice(); //take new AliRun object from galice.root
17 gHL=gAL->GetDetectorLoader("HMPID"); gH=(AliHMPID*)gAL->GetAliRun()->GetDetector("HMPID"); //get HMPID object from galice.root
18 gMaxEvt=gAL->GetNumberOfEvents()-1;
d73ece6a 19 gHL->LoadHits(); gHL->LoadSDigits(); gHL->LoadDigits(); gHL->LoadRecPoints();
8282f9d0 20
21 TFile *pEsdFl=TFile::Open("AliESDs.root"); gEsdTr=(TTree*) pEsdFl->Get("esdTree"); gEsdTr->SetBranchAddress("ESD", &gEsd);
22 pAll->cd(7); TButton *pBtn=new TButton("Next","ReadEvt()",0,0,0.2,0.1); pBtn->Draw();
d73ece6a 23 TButton *pHitBtn=new TButton("Print hits","PrintHits()",0,0.2,0.3,0.3); pHitBtn->Draw();
24 TButton *pSdiBtn=new TButton("Print sdis","PrintSdis()",0,0.4,0.3,0.5); pSdiBtn->Draw();
25 TButton *pDigBtn=new TButton("Print digs","PrintDigs()",0,0.6,0.3,0.7); pDigBtn->Draw();
26 TButton *pCluBtn=new TButton("Print clus","PrintClus()",0,0.8,0.3,0.9); pCluBtn->Draw();
8282f9d0 27 ReadEvt();
28 }else{
29 pAll->cd(7); TButton *pBtn=new TButton("Next","SimEvt()",0,0,0.2,0.1); pBtn->Draw();
30 SimEvt();
31 }
32}
33//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
34void ReadEvt()
35{// Read curent event and display it assumes that session is alredy opened
36 TClonesArray hits("AliHMPIDHit");
37 if(gEvt>gMaxEvt) gEvt=0; if(gEvt<0) gEvt=gMaxEvt;
38
39 gEsdTr->GetEntry(gEvt); gAL->GetEvent(gEvt);
d73ece6a 40 ReadHits(&hits); gHL->TreeS()->GetEntry(0); gHL->TreeD()->GetEntry(0); gHL->TreeR()->GetEntry(0);
41
8282f9d0 42
8282f9d0 43 DrawEvt(&hits,gH->DigLst(),gH->CluLst(),gEsd);
44 gEvt++;
45}
46//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
47void SimEvt()
48{
49 TClonesArray hits("AliHMPIDHit");
50 TClonesArray sdig("AliHMPIDDigit");
51 TObjArray digs(7); for(Int_t i=0;i<7;i++) digs.AddAt(new TClonesArray("AliHMPIDDigit"),i);
52 TObjArray clus(7); for(Int_t i=0;i<7;i++) clus.AddAt(new TClonesArray("AliHMPIDCluster"),i);
53 AliESD esd;
27311693 54 AliHMPIDDigit::fSigmas=4;
55 AliHMPIDDigitizer::DoNoise(kFALSE);
8282f9d0 56 SimEsd(&esd);
57 SimHits(&esd,&hits);
58 AliHMPIDv1::Hit2Sdi(&hits,&sdig);
59 AliHMPIDDigitizer::Sdi2Dig(&sdig,&digs);
50e69536 60 AliHMPIDReconstructor::Dig2Clu(&digs,&clus);
bd73b83f 61 AliHMPIDTracker::Recon(&esd,&clus,pNmean);
8282f9d0 62
63 pAll->cd(3); gPad->Clear(); TLatex txt;txt.DrawLatex(0.2,0.2,Form("Simulated event %i",gEvt));
64 DrawEvt(&hits,&digs,&clus,&esd);
65 gEvt++;
66}//SimEvt()
a1d55ff3 67//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
68void SimEsd(AliESD *pEsd)
69{
70 TParticle part; TLorentzVector mom;
50e69536 71 for(Int_t iTrk=0;iTrk<100;iTrk++){//stack loop
a1d55ff3 72 part.SetPdgCode(kProton);
50e69536 73 part.SetProductionVertex(0,0,0,0);
74 Double_t eta= -0.2+gRandom->Rndm()*0.4; //rapidity is random [-0.2,+0.2]
75 Double_t phi= gRandom->Rndm()*60.*TMath::DegToRad(); //phi is random [ 0 , 60 ] degrees
76 mom.SetPtEtaPhiM(5,eta,phi,part.GetMass());
a1d55ff3 77 part.SetMomentum(mom);
78 AliESDtrack trk(&part);
79 pEsd->AddTrack(&trk);
80 }//stack loop
81}//EsdFromStack()
82//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8282f9d0 83void SimHits(AliESD *pEsd, TClonesArray *pHits)
a1d55ff3 84{
50e69536 85 const Int_t kCerenkov=50000050;
86 const Int_t kFeedback=50000051;
f3bae3e2 87
a1d55ff3 88 AliHMPIDRecon rec;
f3bae3e2 89
a591e55f 90 Int_t hc=0;
a1d55ff3 91 for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){//tracks loop
92 AliESDtrack *pTrk=pEsd->GetTrack(iTrk);
93 Float_t xRa,yRa;
94 Int_t ch=AliHMPIDTracker::IntTrkCha(pTrk,xRa,yRa);
95 if(ch<0) continue; //this track does not hit HMPID
50e69536 96 Float_t beta = pTrk->GetP()/(TMath::Sqrt(pTrk->GetP()*pTrk->GetP()+0.938*0.938));
97 Float_t ckov=TMath::ACos(1./(beta*1.292));
a1d55ff3 98
50e69536 99 Float_t theta,phi,xPc,yPc,;
100 pTrk->GetHMPIDtrk(xPc,yPc,theta,phi);
101 rec.SetTrack(xRa,yRa,theta,phi);
a1d55ff3 102
8282f9d0 103 if(!AliHMPIDDigit::IsInDead(xPc,yPc)) new((*pHits)[hc++]) AliHMPIDHit(ch,200e-9,kProton ,iTrk,xPc,yPc); //mip hit
50e69536 104 Int_t nPhots = (Int_t)(20.*TMath::Power(TMath::Sin(ckov),2)/TMath::Power(TMath::Sin(TMath::ACos(1./1.292)),2));
105 for(int i=0;i<nPhots;i++){
a591e55f 106 TVector2 pos;
107 pos=rec.TracePhot(ckov,gRandom->Rndm()*TMath::TwoPi());
50e69536 108 if(!AliHMPIDDigit::IsInDead(pos.X(),pos.Y())) new((*pHits)[hc++]) AliHMPIDHit(ch,7.5e-9,kCerenkov,iTrk,pos.X(),pos.Y());
a1d55ff3 109 } //photon hits
110 }//tracks loop
8282f9d0 111}//SimHits()
a1d55ff3 112//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8282f9d0 113void ReadHits(TClonesArray *pHitLst)
a1d55ff3 114{
8282f9d0 115 pHitLst->Delete(); Int_t cnt=0;
116 for(Int_t iEnt=0;iEnt<gHL->TreeH()->GetEntries();iEnt++){ //TreeH loop
117 gHL->TreeH()->GetEntry(iEnt); //get current entry (prim)
118 for(Int_t iHit=0;iHit<gH->Hits()->GetEntries();iHit++){ //hits loop
119 AliHMPIDHit *pHit = (AliHMPIDHit*)gH->Hits()->At(iHit); //get current hit
120 new((*pHitLst)[cnt++]) AliHMPIDHit(*pHit);
121 }//hits loop for this entry
122 }//tree entries loop
123
124}//ReadHits()
a1d55ff3 125//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8282f9d0 126void DrawCh(Int_t iCh)
a1d55ff3 127{
128 gPad->Range(-10,-10,AliHMPIDDigit::SizeAllX()+5,AliHMPIDDigit::SizeAllY()+5);
f3bae3e2 129 if(iCh>=0){TLatex txt; txt.SetTextSize(0.1); txt.DrawLatex(-5,-5,Form("%i",iCh));}
a1d55ff3 130
131 for(Int_t iPc=AliHMPIDDigit::kMinPc;iPc<=AliHMPIDDigit::kMaxPc;iPc++){
132 TBox *pBox=new TBox(AliHMPIDDigit::fMinPcX[iPc],AliHMPIDDigit::fMinPcY[iPc],
133 AliHMPIDDigit::fMaxPcX[iPc],AliHMPIDDigit::fMaxPcY[iPc]);
134
8282f9d0 135 pBox->SetFillStyle(0);
a1d55ff3 136 pBox->Draw();
137
a1d55ff3 138 }//PC loop
8282f9d0 139}//DrawCh()
a1d55ff3 140//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8282f9d0 141void DrawEvt(TClonesArray *pHitLst,TObjArray *pDigLst,TObjArray *pCluLst,AliESD *pEsd)
a1d55ff3 142{//draws all the objects of current event in given canvas
a1d55ff3 143 AliHMPIDRecon rec;
d73ece6a 144 TPolyMarker *pTxC[7], *pRin[7]; TMarker *pMip,*pCko,*pFee,*pDig,*pClu;
145 pMip=new TMarker; pMip->SetMarkerColor(kRed); pMip->SetMarkerStyle(kOpenTriangleUp);
146 pCko=new TMarker; pCko->SetMarkerColor(kRed); pCko->SetMarkerStyle(kOpenCircle);
147 pFee=new TMarker; pFee->SetMarkerColor(kRed); pFee->SetMarkerStyle(kOpenDiamond);
148 pDig=new TMarker; pDig->SetMarkerColor(kGreen);pDig->SetMarkerStyle(kOpenSquare);
149 pClu=new TMarker; pClu->SetMarkerColor(kBlue); pClu->SetMarkerStyle(kStar);
a1d55ff3 150 for(Int_t ch=0;ch<7;ch++){
151 pTxC[ch]=new TPolyMarker; pTxC[ch]->SetMarkerStyle(2); pTxC[ch]->SetMarkerColor(kRed); pTxC[ch]->SetMarkerSize(3);
152 pRin[ch]=new TPolyMarker; pRin[ch]->SetMarkerStyle(6); pRin[ch]->SetMarkerColor(kMagenta);
153 }
154
a1d55ff3 155 for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){//tracks loop to collect cerenkov rings and intersection points
156 AliESDtrack *pTrk=pEsd->GetTrack(iTrk);
157 Int_t ch=pTrk->GetHMPIDcluIdx();
158 if(ch<0) continue; //this track does not hit HMPID
159 ch/=1000000;
a591e55f 160 Float_t th,ph,xRad,yRad; pTrk->GetHMPIDtrk(xRad,yRad,th,ph);//get info on current track
161// pTxC[ch]->SetNextPoint(xPc,yPc); //add new intersection point TEMPORARLY DISABLED...no more available in ESD!
a1d55ff3 162
163 Float_t ckov=pTrk->GetHMPIDsignal(); Float_t err=TMath::Sqrt(pTrk->GetHMPIDchi2());
a1d55ff3 164 if(ckov>0){
394d35c1 165 Printf("theta %f phi %f ckov %f",th*TMath::RadToDeg(),ph*TMath::RadToDeg(),ckov);
b4ad85e9 166 rec.SetTrack(xRad,yRad,th,ph);
a591e55f 167 for(Int_t j=0;j<100;j++){
168 TVector2 pos;
169 pos=rec.TracePhot(ckov,j*0.0628);
50e69536 170 if(!AliHMPIDDigit::IsInDead(pos.X(),pos.Y())) pRin[ch]->SetNextPoint(pos.X(),pos.Y());
171 }
a1d55ff3 172 }
173 }//tracks loop
174
b08ebbb0 175 Int_t totHit=pHitLst->GetEntriesFast(),totDig=0,totClu=0,totTxC=0;
176 Int_t totCkov=0, totFeed=0,totMip=0;
a1d55ff3 177 for(Int_t iCh=0;iCh<7;iCh++){//chambers loop
d73ece6a 178 totTxC+=pTxC[iCh]->GetN();
179 totDig+=((TClonesArray*)pDigLst->At(iCh))->GetEntriesFast();
180 totClu+=((TClonesArray*)pCluLst->At(iCh))->GetEntriesFast();
181
a1d55ff3 182 switch(iCh){
8282f9d0 183 case 6: pAll->cd(1); break; case 5: pAll->cd(2); break;
184 case 4: pAll->cd(4); break; case 3: pAll->cd(5); break; case 2: pAll->cd(6); break;
185 case 1: pAll->cd(8); break; case 0: pAll->cd(9); break;
a1d55ff3 186 }
187 gPad->SetEditable(kTRUE); gPad->Clear();
8282f9d0 188 DrawCh(iCh);
d73ece6a 189
8c1bbb8a 190 ((TClonesArray*)pDigLst->At(iCh))->Draw(); //draw digits
191
b08ebbb0 192 totCkov=0;totFeed=0;totMip=0;
8c1bbb8a 193 for(Int_t iHit=0;iHit<pHitLst->GetEntriesFast();iHit++){ // Draw hits
a1d55ff3 194 AliHMPIDHit *pHit=(AliHMPIDHit*)pHitLst->At(iHit);
d73ece6a 195 if(pHit->Ch()==iCh) pHit->Draw();
b08ebbb0 196 if(pHit->Pid()==50000050) totCkov++;
197 else if(pHit->Pid()==50000051) totFeed++;
198 else totMip++;
d73ece6a 199 }
200
8c1bbb8a 201 ((TClonesArray*)pCluLst->At(iCh))->Draw(); //draw clusters
202 pTxC[iCh]->Draw(); //draw intersections
203 pRin[iCh]->Draw("CLP"); //draw rings
a1d55ff3 204 gPad->SetEditable(kFALSE);
205 }//chambers loop
d73ece6a 206
207 pAll->cd(3); gPad->Clear(); TLegend *pLeg=new TLegend(0.2,0.2,0.8,0.8);
208 pLeg->SetHeader(Form("Event %i Total %i",gEvt,gMaxEvt+1));
209 pLeg->AddEntry(pTxC[0],Form("TRKxPC %i" ,totTxC),"p");
b08ebbb0 210 pLeg->AddEntry(pMip ,Form("Mip hits %i" ,totMip),"p");
211 pLeg->AddEntry(pCko ,Form("Ckov hits %i" ,totCkov),"p");
212 pLeg->AddEntry(pFee ,Form("Feed hits %i" ,totFeed),"p");
d73ece6a 213 pLeg->AddEntry(pDig ,Form("Digs %i" ,totDig),"p");
214 pLeg->AddEntry(pClu ,Form("Clus %i" ,totClu),"p");
215 pLeg->Draw();
a1d55ff3 216}//DrawEvt()
8282f9d0 217//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
218void DoZoom(Int_t evt, Int_t px, Int_t py, TObject *obj)
a1d55ff3 219{
8282f9d0 220 if(evt!=5 && evt!=6) return; //5- zoom in 6-zoom out
221 const Int_t minZoom=64;
222 const Int_t maxZoom=2;
223 static Int_t zoom=minZoom; //zoom level
224 if(evt==5&&zoom==maxZoom) return;
225 if(evt==6&&zoom==minZoom) return;
226
d73ece6a 227 // if(!obj->IsA()->InheritsFrom("TPad")) return; //current object is not pad
50e69536 228 TVirtualPad *pPad=gPad->GetSelectedPad();
8282f9d0 229 if(pPad->GetNumber()==3 || pPad->GetNumber()==7) return; //current pad is wrong
230
231 // Printf("evt=%i (%i,%i) %s",evt,px,py,obj->GetName());
232
233 Float_t x=pPad->AbsPixeltoX(px); Float_t y=pPad->AbsPixeltoY(py);
6ed6a312 234
235 if(evt==5){ zoom=zoom/2; pPad->Range(x-zoom*2,y-zoom*2,x+zoom*2,y+zoom*2);} //zoom in
236 else { zoom=zoom*2; pPad->Range(x-zoom*2,y-zoom*2,x+zoom*2,y+zoom*2);} //zoom out
237 if(zoom==minZoom) pPad->Range(-10,-10,AliHMPIDDigit::SizeAllX()+5,AliHMPIDDigit::SizeAllY()+5);
8282f9d0 238 ((TCanvas *)gTQSender)->SetTitle(Form("zoom x%i",minZoom/zoom));
239 pPad->Modified();
240 pPad->Update();
a1d55ff3 241}
8282f9d0 242//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
d73ece6a 243void PrintHits()
244{
245//Prints a list of HMPID hits for a given event. Default is event number 0.
246 Printf("List of HMPID hits for event %i",gEvt);
247
248 Int_t iTotHits=0;
249 for(Int_t iPrim=0;iPrim<gHL->TreeH()->GetEntries();iPrim++){//prims loop
250 gHL->TreeH()->GetEntry(iPrim);
251 gH->Hits()->Print();
252 iTotHits+=gH->Hits()->GetEntries();
253 }
254 Printf("totally %i hits for event %i",iTotHits,gEvt);
255}
256//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
257void PrintSdis()
258{
259//prints a list of HMPID sdigits for a given event
260 Printf("List of HMPID sdigits for event %i",gEvt);
261 gH->SdiLst()->Print();
262 Printf("totally %i sdigits for event %i",gH->SdiLst()->GetEntries(),gEvt);
263}
264//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
265void PrintDigs()
266{
267//prints a list of HMPID digits
268 Printf("List of HMPID digits for event %i",gEvt);
269 gH->DigLst()->Print();
270 Int_t totDigs=0; for(Int_t i=0;i<7;i++) {totDigs+=gH->DigLst(i)->GetEntries();}
271 Printf("totally %i digits for event %i",totDigs,gEvt);
272}
273//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
274void PrintClus()
275{//prints a list of HMPID clusters for a given event
276 Printf("List of HMPID clusters for event %i",gEvt);
277 gH->CluLst()->Print();
278
279 Int_t iCluCnt=0; for(Int_t iCh=0;iCh<7;iCh++) iCluCnt+=gH->CluLst(iCh)->GetEntries();
280
281 Printf("totally %i clusters for event %i",iCluCnt,gEvt);
282}
50e69536 283//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
bd73b83f 284void OpenCalib()
50e69536 285{
286 AliCDBManager* pCDB = AliCDBManager::Instance();
287 pCDB->SetDefaultStorage("local://$HOME");
288 AliCDBEntry *pQthreEnt=pCDB->Get("HMPID/Calib/Qthre",0);
289 AliCDBEntry *pNmeanEnt=pCDB->Get("HMPID/Calib/Nmean",0);
290
291 if(!pQthreEnt || ! pNmeanEnt) return;
292
bd73b83f 293 pNmean=(TObjArray*)pNmeanEnt->GetObject();
50e69536 294}