Theta Cerenkov fixed + minors
[u/mrichter/AliRoot.git] / HMPID / Hdisp.C
1 TCanvas *pAll=0;
2 AliRunLoader *gAL=0; AliLoader *gHL=0; AliESD *gEsd=0; TTree *gEsdTr=0; AliHMPID *gH=0;
3 Int_t gEvt=0; Int_t gMaxEvt=0;
4 TObjArray *pNmean;
5 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6 void Hdisp()
7 {//display events from files if any in current directory or simulated events
8   pAll=new TCanvas("all","",1300,900); pAll->Divide(3,3,0,0);
9 //  pAll->ToggleEditor();
10   pAll->Connect("ProcessedEvent(Int_t,Int_t,Int_t,TObject*)",0,"","DoZoom(Int_t,Int_t,Int_t,TObject*)");
11
12   OpenCalib();  
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;
19     gHL->LoadHits(); gHL->LoadSDigits(); gHL->LoadDigits(); gHL->LoadRecPoints();
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();
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();  
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 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
34 void 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); 
40   ReadHits(&hits); gHL->TreeS()->GetEntry(0); gHL->TreeD()->GetEntry(0); gHL->TreeR()->GetEntry(0);
41   
42   
43   DrawEvt(&hits,gH->DigLst(),gH->CluLst(),gEsd);
44   gEvt++;
45 }
46 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
47 void 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;
54   AliHMPIDDigit::fSigmas=4;
55   AliHMPIDDigitizer::DoNoise(kFALSE);
56   SimEsd(&esd);
57   SimHits(&esd,&hits);
58              AliHMPIDv1::Hit2Sdi(&hits,&sdig);                               
59       AliHMPIDDigitizer::Sdi2Dig(&sdig,&digs);     
60       AliHMPIDReconstructor::Dig2Clu(&digs,&clus);
61       AliHMPIDTracker::Recon(&esd,&clus,pNmean);
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()
67 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
68 void SimEsd(AliESD *pEsd)
69 {
70   TParticle part; TLorentzVector mom;
71   for(Int_t iTrk=0;iTrk<100;iTrk++){//stack loop
72     part.SetPdgCode(kProton);
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());
77     part.SetMomentum(mom);
78     AliESDtrack trk(&part);
79     pEsd->AddTrack(&trk);
80   }//stack loop  
81 }//EsdFromStack()
82 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
83 void SimHits(AliESD *pEsd, TClonesArray *pHits)
84 {
85   const Int_t kCerenkov=50000050;
86   const Int_t kFeedback=50000051;
87   
88   AliHMPIDRecon rec;
89   
90   Int_t hc=0; TVector2 pos;
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
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));
98
99     Float_t theta,phi,xPc,yPc,; 
100     pTrk->GetHMPIDtrk(xPc,yPc,theta,phi); 
101     rec.SetTrack(xRa,yRa,theta,phi); 
102     
103     if(!AliHMPIDDigit::IsInDead(xPc,yPc)) new((*pHits)[hc++]) AliHMPIDHit(ch,200e-9,kProton  ,iTrk,xPc,yPc);                 //mip hit
104     /*
105     for(int i=0;i<4;i++) { 
106       Float_t x=gRandom->Rndm()*130;Float_t y=gRandom->Rndm()*126;
107       if(!AliHMPIDDigit::IsInDead(x,y)) new((*pHits)[hc++]) AliHMPIDHit(ch,7.5e-9,kFeedback,iTrk,x,y); //bkg hits 4 per track
108     }
109     */
110     Int_t nPhots = (Int_t)(20.*TMath::Power(TMath::Sin(ckov),2)/TMath::Power(TMath::Sin(TMath::ACos(1./1.292)),2));
111     for(int i=0;i<nPhots;i++){
112       rec.TracePhot(ckov,gRandom->Rndm()*TMath::TwoPi(),pos);
113       if(!AliHMPIDDigit::IsInDead(pos.X(),pos.Y())) new((*pHits)[hc++]) AliHMPIDHit(ch,7.5e-9,kCerenkov,iTrk,pos.X(),pos.Y());
114     }                      //photon hits  
115   }//tracks loop    
116 }//SimHits()
117 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
118 void ReadHits(TClonesArray *pHitLst)
119 {
120   pHitLst->Delete();  Int_t cnt=0;
121   for(Int_t iEnt=0;iEnt<gHL->TreeH()->GetEntries();iEnt++){       //TreeH loop
122     gHL->TreeH()->GetEntry(iEnt);                                 //get current entry (prim)                
123     for(Int_t iHit=0;iHit<gH->Hits()->GetEntries();iHit++){       //hits loop
124       AliHMPIDHit *pHit = (AliHMPIDHit*)gH->Hits()->At(iHit);     //get current hit        
125       new((*pHitLst)[cnt++]) AliHMPIDHit(*pHit);
126     }//hits loop for this entry
127   }//tree entries loop
128   
129 }//ReadHits()
130 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
131 void DrawCh(Int_t iCh) 
132
133   gPad->Range(-10,-10,AliHMPIDDigit::SizeAllX()+5,AliHMPIDDigit::SizeAllY()+5); 
134   if(iCh>=0){TLatex txt; txt.SetTextSize(0.1); txt.DrawLatex(-5,-5,Form("%i",iCh));}
135   
136   for(Int_t iPc=AliHMPIDDigit::kMinPc;iPc<=AliHMPIDDigit::kMaxPc;iPc++){
137     TBox *pBox=new TBox(AliHMPIDDigit::fMinPcX[iPc],AliHMPIDDigit::fMinPcY[iPc],
138                         AliHMPIDDigit::fMaxPcX[iPc],AliHMPIDDigit::fMaxPcY[iPc]);
139     
140     pBox->SetFillStyle(0);
141     pBox->Draw();
142     
143   }//PC loop      
144 }//DrawCh()
145 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
146 void DrawEvt(TClonesArray *pHitLst,TObjArray *pDigLst,TObjArray *pCluLst,AliESD *pEsd)
147 {//draws all the objects of current event in given canvas
148   AliHMPIDRecon rec;  
149   TPolyMarker *pTxC[7], *pRin[7]; TMarker *pMip,*pCko,*pFee,*pDig,*pClu;
150   pMip=new TMarker; pMip->SetMarkerColor(kRed);  pMip->SetMarkerStyle(kOpenTriangleUp);
151   pCko=new TMarker; pCko->SetMarkerColor(kRed);  pCko->SetMarkerStyle(kOpenCircle);
152   pFee=new TMarker; pFee->SetMarkerColor(kRed);  pFee->SetMarkerStyle(kOpenDiamond);
153   pDig=new TMarker; pDig->SetMarkerColor(kGreen);pDig->SetMarkerStyle(kOpenSquare);
154   pClu=new TMarker; pClu->SetMarkerColor(kBlue); pClu->SetMarkerStyle(kStar);
155   for(Int_t ch=0;ch<7;ch++){
156     pTxC[ch]=new TPolyMarker; pTxC[ch]->SetMarkerStyle(2); pTxC[ch]->SetMarkerColor(kRed); pTxC[ch]->SetMarkerSize(3);
157     pRin[ch]=new TPolyMarker; pRin[ch]->SetMarkerStyle(6); pRin[ch]->SetMarkerColor(kMagenta);
158   }
159   
160   for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){//tracks loop to collect cerenkov rings and intersection points
161     AliESDtrack *pTrk=pEsd->GetTrack(iTrk);
162     Int_t ch=pTrk->GetHMPIDcluIdx();
163     if(ch<0) continue; //this track does not hit HMPID
164     ch/=1000000; 
165     Float_t th,ph,xPc,yPc; pTrk->GetHMPIDtrk(xPc,yPc,th,ph);  //get info on current track
166     pTxC[ch]->SetNextPoint(xPc,yPc);                          //add new intersection point
167     
168     Float_t ckov=pTrk->GetHMPIDsignal();  Float_t err=TMath::Sqrt(pTrk->GetHMPIDchi2());
169     if(ckov>0){
170       Printf("theta %f phi %f ckov %f",th*TMath::RadToDeg(),ph*TMath::RadToDeg(),ckov);
171       rec.SetTrack(xPc,yPc,th,ph+TMath::Pi());
172       TVector2 pos;
173        Double_t allGapz=rec.fgkWinThick+0.5*rec.fgkRadThick+rec.fgkGapThick;   //to semplify, the (x,y) are calculted from (xPc,yPc) back to radiator in straight line (Bz=0)
174       TVector3 dir(0,0,-1);TVector3 miprad(xPc,yPc,allGapz);
175       AliHMPIDRecon::Propagate(dir,miprad,0);
176       Double_t xRad=miprad.X();
177       Double_t yRad=miprad.Y();
178       for(int j=0;j<100;j++){ 
179         rec.TracePhot(xRad,yRad,ckov,j*0.0628,pos);
180        if(!AliHMPIDDigit::IsInDead(pos.X(),pos.Y())) pRin[ch]->SetNextPoint(pos.X(),pos.Y());
181       }      
182     }
183   }//tracks loop
184       
185   Int_t totHit=pHitLst->GetEntriesFast(),totDig=0,totClu=0,totTxC=0;
186   Int_t totCkov=0, totFeed=0,totMip=0;
187   for(Int_t iCh=0;iCh<7;iCh++){//chambers loop    
188     totTxC+=pTxC[iCh]->GetN();
189     totDig+=((TClonesArray*)pDigLst->At(iCh))->GetEntriesFast();
190     totClu+=((TClonesArray*)pCluLst->At(iCh))->GetEntriesFast();
191     
192     switch(iCh){
193       case 6: pAll->cd(1); break; case 5: pAll->cd(2); break;
194       case 4: pAll->cd(4); break; case 3: pAll->cd(5); break; case 2: pAll->cd(6); break;
195                                   case 1: pAll->cd(8); break; case 0: pAll->cd(9); break;
196     }
197     gPad->SetEditable(kTRUE); gPad->Clear(); 
198     DrawCh(iCh);
199                            
200     ((TClonesArray*)pDigLst->At(iCh))->Draw();               //draw digits
201     
202     totCkov=0;totFeed=0;totMip=0;
203     for(Int_t iHit=0;iHit<pHitLst->GetEntriesFast();iHit++){ // Draw hits
204       AliHMPIDHit *pHit=(AliHMPIDHit*)pHitLst->At(iHit);
205       if(pHit->Ch()==iCh) pHit->Draw();
206       if(pHit->Pid()==50000050) totCkov++;
207       else if(pHit->Pid()==50000051) totFeed++;
208       else totMip++;
209     }    
210            
211     ((TClonesArray*)pCluLst->At(iCh))->Draw();              //draw clusters
212                             pTxC[iCh]->Draw();              //draw intersections
213                             pRin[iCh]->Draw("CLP");         //draw rings
214     gPad->SetEditable(kFALSE);
215   }//chambers loop
216   
217   pAll->cd(3);  gPad->Clear(); TLegend *pLeg=new TLegend(0.2,0.2,0.8,0.8);
218                                         pLeg->SetHeader(Form("Event %i Total %i",gEvt,gMaxEvt+1));
219                                         pLeg->AddEntry(pTxC[0],Form("TRKxPC %i"     ,totTxC),"p");
220                                         pLeg->AddEntry(pMip   ,Form("Mip hits %i"   ,totMip),"p");    
221                                         pLeg->AddEntry(pCko   ,Form("Ckov hits %i"  ,totCkov),"p");    
222                                         pLeg->AddEntry(pFee   ,Form("Feed hits %i"  ,totFeed),"p");    
223                                         pLeg->AddEntry(pDig   ,Form("Digs %i"       ,totDig),"p");    
224                                         pLeg->AddEntry(pClu   ,Form("Clus %i"       ,totClu),"p");    
225                                         pLeg->Draw();
226 }//DrawEvt()
227 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
228 void DoZoom(Int_t evt, Int_t px, Int_t py, TObject *obj)
229 {
230   if(evt!=5 && evt!=6) return; //5- zoom in 6-zoom out
231   const Int_t minZoom=64;
232   const Int_t maxZoom=2;
233   static Int_t zoom=minZoom; //zoom level
234   if(evt==5&&zoom==maxZoom) return; 
235   if(evt==6&&zoom==minZoom) return; 
236   
237  // if(!obj->IsA()->InheritsFrom("TPad")) return;  //current object is not pad
238   TVirtualPad *pPad=gPad->GetSelectedPad();
239   if(pPad->GetNumber()==3 || pPad->GetNumber()==7) return; //current pad is wrong
240
241  // Printf("evt=%i (%i,%i) %s",evt,px,py,obj->GetName());
242     
243   Float_t x=pPad->AbsPixeltoX(px); Float_t y=pPad->AbsPixeltoY(py); 
244  
245   if(evt==5){ zoom=zoom/2;     pPad->Range(x-zoom*2,y-zoom*2,x+zoom*2,y+zoom*2);} //zoom in
246   else      { zoom=zoom*2;     pPad->Range(x-zoom*2,y-zoom*2,x+zoom*2,y+zoom*2);} //zoom out 
247   if(zoom==minZoom) pPad->Range(-10,-10,AliHMPIDDigit::SizeAllX()+5,AliHMPIDDigit::SizeAllY()+5);
248   ((TCanvas *)gTQSender)->SetTitle(Form("zoom x%i",minZoom/zoom));
249   pPad->Modified();
250   pPad->Update();                                              
251 }
252 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
253 void PrintHits()
254 {
255 //Prints a list of HMPID hits for a given event. Default is event number 0.
256   Printf("List of HMPID hits for event %i",gEvt);
257   
258   Int_t iTotHits=0;
259   for(Int_t iPrim=0;iPrim<gHL->TreeH()->GetEntries();iPrim++){//prims loop
260     gHL->TreeH()->GetEntry(iPrim);      
261     gH->Hits()->Print();
262     iTotHits+=gH->Hits()->GetEntries();
263   }
264   Printf("totally %i hits for event %i",iTotHits,gEvt);
265 }
266 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
267 void PrintSdis()
268 {
269 //prints a list of HMPID sdigits  for a given event
270   Printf("List of HMPID sdigits for event %i",gEvt);
271   gH->SdiLst()->Print();
272   Printf("totally %i sdigits for event %i",gH->SdiLst()->GetEntries(),gEvt);
273 }
274 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
275 void PrintDigs()
276 {
277 //prints a list of HMPID digits
278   Printf("List of HMPID digits for event %i",gEvt);  
279   gH->DigLst()->Print();
280   Int_t totDigs=0;  for(Int_t i=0;i<7;i++) {totDigs+=gH->DigLst(i)->GetEntries();}
281   Printf("totally %i digits for event %i",totDigs,gEvt);
282 }
283 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
284 void PrintClus()
285 {//prints a list of HMPID clusters  for a given event
286   Printf("List of HMPID clusters for event %i",gEvt);
287   gH->CluLst()->Print();
288   
289   Int_t iCluCnt=0; for(Int_t iCh=0;iCh<7;iCh++) iCluCnt+=gH->CluLst(iCh)->GetEntries();
290   
291   Printf("totally %i clusters for event %i",iCluCnt,gEvt);
292 }
293 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
294 void OpenCalib()
295 {
296   AliCDBManager* pCDB = AliCDBManager::Instance();
297   pCDB->SetDefaultStorage("local://$HOME");
298   AliCDBEntry *pQthreEnt=pCDB->Get("HMPID/Calib/Qthre",0);
299   AliCDBEntry *pNmeanEnt=pCDB->Get("HMPID/Calib/Nmean",0);
300   
301   if(!pQthreEnt || ! pNmeanEnt) return;
302   
303   pNmean=(TObjArray*)pNmeanEnt->GetObject(); 
304 }