First prototype of the new AliCluster base class
[u/mrichter/AliRoot.git] / HMPID / Hdisp.C
1 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 void SimEsd(AliESD *pEsd)
3 {
4   TParticle part; TLorentzVector mom;
5   for(Int_t iTrk=0;iTrk<25;iTrk++){//stack loop
6     part.SetPdgCode(kProton);
7     part.SetProductionVertex(0,0,0,0);  
8     Double_t eta= -0.4+gRandom->Rndm()*0.8; //rapidity is random [-0.4,+0.4]
9     Double_t phi= gRandom->Rndm()*1.4;      //phi is random      [ 0  , 80 ] degrees    
10     mom.SetPtEtaPhiM(2,eta,phi,part.GetMass());
11     part.SetMomentum(mom);
12     AliESDtrack trk(&part);
13     pEsd->AddTrack(&trk);
14   }//stack loop  
15 }//EsdFromStack()
16 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
17 void EsdFromStack(AliESD *pEsd)
18 {
19   al->LoadHeader();al->LoadKinematics();
20   AliStack *pStk=al->Stack();  
21   
22   for(Int_t iTrk=0;iTrk<pStk->GetNtrack();iTrk++){//stack loop
23     TParticle *pPart=pStk->Particle(iTrk);
24     if(pPart->GetPDG()->Charge()==0) continue; //neutral particles are not reconstructed
25     if(pPart->GetFirstMother()>0)    continue; //do not consider secondaries
26     AliESDtrack trk(pPart);
27     pEsd->AddTrack(&trk);
28   }//stack loop  
29 }//EsdFromStack()
30 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
31 void HitsFromEsd(AliESD *pEsd, TClonesArray *pHitLst)
32 {
33   AliHMPIDRecon rec;
34   const Int_t kCerenkov=50000050,kFeedback=50000051;
35   Int_t hc=0; TVector2 pos;
36   for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){//tracks loop
37     AliESDtrack *pTrk=pEsd->GetTrack(iTrk);
38     Float_t xRa,yRa;
39     Int_t ch=AliHMPIDTracker::IntTrkCha(pTrk,xRa,yRa);
40     if(ch<0) continue; //this track does not hit HMPID
41     Float_t ckov=0.63;
42
43     Float_t th,ph,xPc,yPc,; pTrk->GetHMPIDtrk(xPc,yPc,th,ph); rec.SetTrack(xRa,yRa,th,ph); 
44     
45     if(!AliHMPIDDigit::IsInDead(xPc,yPc)) new((*pHitLst)[hc++]) AliHMPIDHit(ch,200e-9,kProton  ,iTrk,xPc,yPc);                 //mip hit
46     for(int i=0;i<4;i++)  new((*pHitLst)[hc++]) AliHMPIDHit(ch,7.5e-9,kFeedback,iTrk,gRandom->Rndm()*130,gRandom->Rndm()*126); //bkg hits 4 per track
47     for(int i=0;i<16;i++){
48       rec.TracePhot(ckov,gRandom->Rndm()*TMath::TwoPi(),pos);
49       new((*pHitLst)[hc++]) AliHMPIDHit(ch,7.5e-9,kCerenkov,iTrk,pos.X(),pos.Y());
50     }                      //photon hits  
51   }//tracks loop    
52 }//HitsFromEsd()
53 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
54 void DrawSeg()
55 {
56   TCanvas *pC=new TCanvas("seg","Segmentation as seen from electronics side");
57   DrawPc(1);
58
59   pC->ToggleEventStatus();
60   pC->SetEditable(0);
61   pC->AddExec("status","DrawStatus()");
62 }
63 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
64 void DrawStatus()
65 {// Show info about current cursur position in status bar of the canvas
66 //  Printf("event %i",gPad->GetEvent()); return;
67   TPad *pad=gPad;
68   TCanvas *pC=(TCanvas*)pad; 
69   TRootCanvas *pRC= (TRootCanvas*)pC->GetCanvasImp();
70   TGStatusBar *pBar=pRC->GetStatusBar();
71   pBar->SetParts(5);
72   Float_t x=pad->AbsPixeltoX(pad->GetEventX()); Float_t y=pad->AbsPixeltoY(pad->GetEventY());
73   if(AliHMPIDDigit::IsInDead(x,y))
74     pBar->SetText("Out of sensitive area",4);    
75   else{
76     Int_t pc,px,py,w32,ddl,r,d,a;  AliHMPIDDigit::Lors2Pad(x,y,pc,px,py); AliHMPIDDigit dig; dig.Set(1,pc,px,py); dig.Raw(w32,ddl,r,d,a);
77     pBar->SetText(Form("(pc%i,px%i,py%i) (r%i,d%i,a%i) (%.2f,%.2f)",
78                          pc  ,px  ,py,    r  ,d  ,a   ,dig.LorsX(),dig.LorsY()),4);
79   }    
80 //  if(pad->GetEvent()==1){
81 //    new TCanvas("zoom",Form("Row %i DILOGIC %i",dig.Row(),dig.Dilogic()));  
82 //  }
83 }//DrawStatus()
84 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
85 void DrawPc(Bool_t isFill) 
86
87   gPad->Range(-10,-10,AliHMPIDDigit::SizeAllX()+5,AliHMPIDDigit::SizeAllY()+5); 
88   
89     
90   Float_t dX2=0.5*AliHMPIDDigit::SizePadX(),
91           dY2=0.5*AliHMPIDDigit::SizePadY() ;
92   
93   TLatex txt; txt.SetTextSize(0.01);
94   TLine *pL;
95   
96   AliHMPIDDigit dig;   UInt_t w32; Int_t ddl,r,d,a;    
97   
98   for(Int_t iPc=AliHMPIDDigit::kMinPc;iPc<=AliHMPIDDigit::kMaxPc;iPc++){
99     TBox *pBox=new TBox(AliHMPIDDigit::fMinPcX[iPc],AliHMPIDDigit::fMinPcY[iPc],
100                         AliHMPIDDigit::fMaxPcX[iPc],AliHMPIDDigit::fMaxPcY[iPc]);
101     
102     if(iPc==0||iPc==2||iPc==4) pBox->SetFillColor(29);
103     else                       pBox->SetFillColor(41);
104     pBox->Draw();
105     
106     if(!isFill)  continue;
107     
108 //    if(iPc%2) {dig.Set(0,iPc,79,25); txt.DrawText(dig.LorsX()+2,dig.LorsY(),Form("PC%i",dig.Pc()));}//print PC#    
109
110     txt.SetTextAlign(32);
111     for(Int_t iRow=0;iRow<8 ;iRow++){//draw row lines (horizontal)
112       dig.Set(0,iPc,0,iRow*6); dig.Raw(w32,ddl,r,d,a);  //set digit to the left-down pad of this row
113       
114       if(iPc%2) txt.DrawText(dig.LorsX()-1           ,dig.LorsY(),Form("%i",dig.PadPcY()));                  //print PadY#    
115                 txt.DrawText(dig.LorsX()-1+(iPc%2)*67,dig.LorsY()+2,Form("r%i",r));                          //print Row#    
116       pL=new TLine(dig.LorsX()-dX2,dig.LorsY()-dY2,dig.LorsX()+AliHMPIDDigit::SizePcX()-dX2,dig.LorsY()-dY2);//draw horizontal line 
117       if(iRow!=0) pL->Draw(); 
118     }//row loop  
119     
120     txt.SetTextAlign(13);
121     for(Int_t iDil=0;iDil<10;iDil++){//draw dilogic lines (vertical)
122       dig.Set(0,iPc,iDil*8,0); dig.Raw(w32,ddl,r,d,a);       //set this digit to the left-down pad of this dilogic        
123       
124                            txt.DrawText(dig.LorsX()  ,dig.LorsY()-1,Form("%i",dig.PadPcX()));                 //print PadX# 
125       if(iPc==4 || iPc==5) txt.DrawText(dig.LorsX()+2,dig.LorsY()+42,Form("d%i",d));              //print Dilogic#    
126       pL=new TLine(dig.LorsX()-dX2,dig.LorsY()-dY2,dig.LorsX()-dX2,dig.LorsY()+AliHMPIDDigit::SizePcY()-dY2); //draw vertical line
127       if(iDil!=0)pL->Draw();
128     }//dilogic loop        
129   }//PC loop      
130 }//DrawPc()
131 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
132 void Close_hed()
133 {
134   TCanvas *pC = ((TCanvas*)gROOT->FindObject("hed"));if(!pC) return;
135   pC->Close();
136   pC=0x0;
137 }
138 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
139 void ReadHits(TClonesArray *pHitLst)
140 {
141   pHitLst->Delete();  Int_t cnt=0;
142   for(Int_t iEnt=0;iEnt<hl->TreeH()->GetEntries();iEnt++){       //TreeH loop
143     hl->TreeH()->GetEntry(iEnt);                                 //get current entry (prim)                
144     for(Int_t iHit=0;iHit<h->Hits()->GetEntries();iHit++){       //hits loop
145       AliHMPIDHit *pHit = (AliHMPIDHit*)h->Hits()->At(iHit);     //get current hit        
146       new((*pHitLst)[cnt++]) AliHMPIDHit(*pHit);
147     }//hits loop for this entry
148   }//tree entries loop
149   
150 }//ReadHits()
151 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
152 void sed()
153 {
154
155   static TCanvas *pC1=0;
156   
157   if(!pC1){
158     pC1=new TCanvas("hed","Simulated evets-View from electronics side, IP is behind the picture.",1000,900); pC1->Divide(3,3);
159     pC1->cd(7); TButton *pBtn=new TButton("Next","sed()",0,0,0.2,0.1);   pBtn->Draw(); 
160     pC1->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", 0, "","zoom(Int_t,Int_t,Int_t,TObject*)");
161   }
162
163   TClonesArray lh("AliHMPIDHit"); 
164   TClonesArray ls("AliHMPIDDigit"); 
165   TObjArray    ld(7); for(Int_t i=0;i<7;i++) ld.AddAt(new TClonesArray("AliHMPIDDigit"),i);
166   TObjArray    lc(7); for(Int_t i=0;i<7;i++) lc.AddAt(new TClonesArray("AliHMPIDCluster"),i);
167   AliESD esd;
168   
169   
170   SimEsd(&esd);
171   HitsFromEsd(&esd,&lh);
172              AliHMPIDv1::Hit2Sdi(&lh,&ls);                               
173       AliHMPIDDigitizer::Sdi2Dig(&ls,&ld);     
174   AliHMPIDReconstructor::Dig2Clu(&ld,&lc);
175 //        AliHMPIDTracker::Recon(&esd,&cl);
176   
177   DrawEvt(pC1,&lh,&ld,&lc,&esd);  
178 }//SimEvt()
179 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
180 void zoom(Int_t evt, Int_t pixx, Int_t pixy, TObject *obj)
181 {
182   Printf("pad %i",gPad->GetNumber());
183   static Int_t lvl=32; //zoom level
184   if(evt!=5 && evt!=6) return;
185   if(evt==5&&lvl==2)  return; //max zoom in
186   if(evt==6&&lvl==32) return; //max zoom out
187   
188   Float_t x=gPad->AbsPixeltoX(pixx); Float_t y=gPad->AbsPixeltoY(pixy); 
189              
190   if(evt==5){ lvl=lvl/2; gPad->Range(x-lvl*2,y-lvl*2,x+lvl*2,y+lvl*2);} //zoom in
191   else      { lvl=32;    gPad->Range(-10,-10,150,140); } //zoom out 
192   ((TCanvas *)gTQSender)->SetTitle(Form("zoom %i",lvl));
193   gPad->Modified();
194   gPad->Update();                                              
195 }
196 void hed()
197 {//event display from files
198   static TCanvas *pC=0;
199   static Int_t iEvt=0;
200   static Int_t iEvtTot=999;
201   static TFile *pEsdFl=0;
202   static TTree *pEsdTr=0;
203   static AliESD *pEsd=0;
204   
205   
206
207   
208   if(!pC&&iEvt<iEvtTot){
209     iEvt=0;
210     iEvtTot=999;
211     if(hl==0) {Printf("hed: no HMPID loader");return;}
212     Printf("Opening session");
213     pEsdFl=TFile::Open("AliESDs.root");     if(!pEsdFl || !pEsdFl->IsOpen()) return;//open AliESDs.root
214     pEsdTr=(TTree*) pEsdFl->Get("esdTree"); if(!pEsdTr)                      return;//get ESD tree
215     pEsdTr->SetBranchAddress("ESD", &pEsd);
216     hl->LoadHits(); hl->LoadDigits(); hl->LoadRecPoints();
217     iEvtTot=pEsdTr->GetEntries();
218     pC=new TCanvas("hed","View from electronics side, IP is behind the picture.",1000,900);  pC->ToggleEventStatus(); pC->Divide(3,3);
219     pC->cd(7); TButton *pBtn=new TButton("Next","hed()",0,0,0.2,0.1);   pBtn->Draw();
220     pC->cd(7); TButton *pBtn=new TButton("Quit","Close_hed()",0.2,0,0.4,0.1);   pBtn->Draw(); 
221     new TGedEditor(pC);
222   }
223  
224   TLatex txt; txt.SetTextSize(0.1);
225   TClonesArray hits("AliHMPIDHit");
226       
227   if(iEvt<iEvtTot){
228     pEsdTr->GetEntry(iEvt); al->GetEvent(iEvt); 
229     hl->TreeD()->GetEntry(0); hl->TreeR()->GetEntry(0);
230     ReadHits(&hits); 
231      
232     pC->cd(3);  gPad->Clear(); txt.DrawLatex(0.2,0.2,Form("Event %i (total %i)",iEvt,iEvtTot));
233     DrawEvt(pC,&hits,h->DigLst(),h->CluLst(),pEsd);
234     
235     iEvt++;
236   }else{
237     Printf("--- No more events available...Bye.");
238     pC->Close();
239     pC=0x0;
240     iEvt=0;
241     iEvtTot=999;
242   }
243 }//hed()
244
245
246
247
248
249
250
251
252
253
254
255
256
257 void DrawEvt(TCanvas *pC,TClonesArray *pHitLst,TObjArray *pDigLst,TObjArray *pCluLst,AliESD *pEsd)
258 {//draws all the objects of current event in given canvas
259
260   AliHMPIDRecon rec;  
261   TPolyMarker *pTxC[7];  TPolyMarker *pRin[7]; //intesections and rings
262   for(Int_t ch=0;ch<7;ch++){
263     pTxC[ch]=new TPolyMarker; pTxC[ch]->SetMarkerStyle(2); pTxC[ch]->SetMarkerColor(kRed); pTxC[ch]->SetMarkerSize(3);
264     pRin[ch]=new TPolyMarker; pRin[ch]->SetMarkerStyle(6); pRin[ch]->SetMarkerColor(kMagenta);
265   }
266   
267   
268   for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){//tracks loop to collect cerenkov rings and intersection points
269     AliESDtrack *pTrk=pEsd->GetTrack(iTrk);
270     Int_t ch=pTrk->GetHMPIDcluIdx();
271     if(ch<0) continue; //this track does not hit HMPID
272     ch/=1000000; 
273     Float_t th,ph,xPc,yPc; pTrk->GetHMPIDtrk(xPc,yPc,th,ph);  //get info on current track
274     pTxC[ch]->SetNextPoint(xPc,yPc);                          //add new intersection point
275     
276     Float_t ckov=pTrk->GetHMPIDsignal();  Float_t err=TMath::Sqrt(pTrk->GetHMPIDchi2());
277     
278     if(ckov>0){
279       rec.SetTrack(xPc,yPc,th,ph);
280      TVector2 pos;  for(int j=0;j<100;j++){rec.TracePhot(ckov,j*0.0628,pos); pRin[ch]->SetNextPoint(pos.X(),pos.Y());}      
281     }
282   }//tracks loop
283       
284   for(Int_t iCh=0;iCh<7;iCh++){//chambers loop    
285     switch(iCh){
286       case 6: pC->cd(1); break; case 5: pC->cd(2); break;
287       case 4: pC->cd(4); break; case 3: pC->cd(5); break; case 2: pC->cd(6); break;
288                                 case 1: pC->cd(8); break; case 0: pC->cd(9); break;
289     }
290     gPad->SetEditable(kTRUE); gPad->Clear(); 
291     DrawPc(0);
292     for(Int_t iHit=0;iHit<pHitLst->GetEntries();iHit++){
293       AliHMPIDHit *pHit=(AliHMPIDHit*)pHitLst->At(iHit);
294       if(pHit->Ch()==iCh)        pHit->Draw();  //draw hits
295     }
296     ((TClonesArray*)pDigLst->At(iCh))->Draw();  //draw digits
297     ((TClonesArray*)pCluLst->At(iCh))->Draw();  //draw clusters
298                             pTxC[iCh]->Draw();  //draw intersections
299                             pRin[iCh]->Draw();  //draw rings
300     gPad->SetEditable(kFALSE);
301   }//chambers loop
302   
303   
304 //  TLatex txt; txt.SetTextSize(0.02);
305 //  txt.DrawLatex(20,-5,Form("#theta=%.4f#pm%.5f with %2i #check{C}"          ,simCkov,simErr,simN));
306 //  txt.DrawLatex(25,-5,Form("#theta=%.4f#pm%.5f with %2i #check{C}"          ,recCkov,recErr,recN));
307 //  txt.DrawLatex(0 ,127,Form("#theta=%.2f#circ   #phi=%.2f#circ @(%.2f,%.2f) ",th*TMath::RadToDeg(),ph*TMath::RadToDeg(),radx,rady));
308 //  Printf("DIG------DIG---------DIG--------DIG------DIG------DIG");pDigLst->Print();Printf("");                   
309 }//DrawEvt()
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350 Reve()
351 {
352   gSystem->Load("libMinuit.so");
353   gSystem->Load("libVMC.so");
354   gSystem->Load("libESD.so");
355   gSystem->Load("libSTEER.so");
356   gSystem->Load("libCDB.so");
357
358   gSystem->Load("libGed.so");
359   gSystem->Load("libRGL.so");
360   gSystem->Load("libGeom.so");
361   gSystem->Load("libVMC.so");
362
363   gSystem->Load("libReve.so");
364
365   gSystem->Load("libHMPIDbase.so");
366   gSystem->Load("libHMPIDsim.so");
367   gSystem->Load("libHMPIDrec.so");
368   
369   TGeoManager::Import("geometry.root");
370   
371   AliHMPIDParam *pParam=AliHMPIDParam::Instance();
372   
373   AliRunLoader *pAL=AliRunLoader::Open(); pAL->LoadHits  ("HMPID"); TTree *pHitT=pAL->GetTreeH("HMPID", false);
374                                           pAL->LoadDigits("HMPID"); TTree *pDigT=pAL->GetTreeD("HMPID", false); 
375                                           
376   Reve::PointSet *pHitPnt=new Reve::PointSet("Hits"));
377   TPolyMarker3D  *pDigPnt=new TPolyMarker3D;  pDigPnt->SetName("Digits"); pDigPnt->SetMarkerColor(kGreen);iPntCnt=0;
378                                           
379   TPointSelector ps(pHitT,pHitPnt,"fX:fY:fZ",""); ps.Select();
380   
381   TClonesArray *pDigLst=new TClonesArray("AliHMPIDDigit"); //this is tmp dig list per chamber
382  
383   for(Int_t iCh=0;iCh<7;iCh++){
384     pDigT->SetBranchAddress(Form("HMPID%i",iCh+1),&pDigLst);
385     pDigT->GetEntry(0);
386     for(Int_t iDig=0;iDig<pDigLst->GetEntries();iDig++){
387       AliHMPIDDigit *pDig=(AliHMPIDDigit*)pDigLst->At(iDig);    
388       TVector2 lors=pParam->Pad2Loc(pDig->PadX(),pDig->PadY());
389       TVector3 mars=pParam->Lors2Mars(iCh,lors.X(),lors.Y());
390       pDigPnt->SetPoint(iPntCnt++,mars.X(),mars.Y(),mars.Z());
391     }//digits loop for chamber        
392   }//chambers loop
393   
394   if(!gReve) new Reve::RGTopFrame(0,0,0,2);
395   gReve->AddGlobalRenderElement(new Reve::RenderElementObjPtr(pDigPnt));
396   gReve->AddRenderElement(pHitPnt);
397   gReve->Redraw3D();
398 }
399
400
401
402 void zzz(Int_t evt,Int_t px,Int_t py,TObject*o)
403 {
404   Printf("evt %i (%i,%i) class %s",evt,px,py,o->IsA()->ClassName());
405 }
406
407 void Gui()
408 {
409   TGMainFrame   *pMF     =new TGMainFrame(gClient->GetRoot(),300,400);//main frame
410 //1 level widgets: button and 2 horizontal frames  
411   TRootEmbeddedCanvas *pDis;
412   pMF->AddFrame(pDis=new TRootEmbeddedCanvas("display",pMF,800,650));
413   pDis->GetCanvas()->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)", 0, "","zzz(Int_t,Int_t,Int_t,TObject*)");
414   
415   pMF->Layout();
416   pMF->MapSubwindows();
417   pMF->Resize(pMF->GetDefaultSize());
418   pMF->MapWindow();
419 }