Cluster improvements + drawing of digits like TBox.
[u/mrichter/AliRoot.git] / HMPID / Hmenu.C
1 AliRun     *a; AliRunLoader *al;   TGeoManager *g; //globals for easy manual manipulations
2 AliHMPID   *h; AliLoader    *hl; AliHMPIDParam *hp;
3 Bool_t isGeomType=kFALSE;
4
5 Int_t nEvent=0;
6
7 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
8 void GetParam()
9 {
10   isGeomType=!isGeomType;
11   if(g) delete g;  if(hp) delete hp; //delete current TGeoManager and AliHMPIDParam
12   if(isGeomType) g=TGeoManager::Import("geometry.root");
13   else           g=TGeoManager::Import("misaligned_geometry.root");
14   hp=AliHMPIDParam::Instance();
15 }//GetParam()
16 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
17 void Hmenu()
18 {   
19   TString status="Status: ";
20   if(gSystem->IsFileInIncludePath("galice.root")){
21     status+="galice.root found";
22     al=AliRunLoader::Open();                                                //try to open galice.root from current dir 
23     if(gAlice) delete gAlice;                                               //in case we execute this in aliroot delete default AliRun object 
24     al->LoadgAlice(); a=al->GetAliRun();                                    //take new AliRun object from galice.root   
25     hl=al->GetDetectorLoader("HMPID");  h=(AliHMPID*)a->GetDetector("HMPID");  //get HMPID object from galice.root
26     
27     status+=Form(" with %i event(s)",al->GetNumberOfEvents()); status+=(h)? " with HMPID": " without HMPID";
28   }else  
29     status+="No galice.root";
30   
31   GetParam();
32   TControlBar *pMenu = new TControlBar("horizontal",status.Data(),0,0);
33     pMenu->AddButton("                     ","","");
34     pMenu->AddButton("       General       ","General()"  ,"general items which do not depend on any files");
35     pMenu->AddButton("                     ",""           ,"");
36     pMenu->AddButton("       Sim data      ","SimData()"  ,"items which expect to have simulated files"    );
37     pMenu->AddButton("                     ",""           ,"");
38     pMenu->AddButton("       Raw data      ","RawData()"  ,"items which expect to have raw files"          );
39     pMenu->AddButton("                     ","       "    ,"");
40     pMenu->AddButton("         Test        ","Test()"     ,"all test utilities");
41     pMenu->AddButton("      PREV EVENT     ","PrevEvent()" ,"Set the previous event"             );
42     pMenu->AddButton("      NEXT EVENT     ","NextEvent()","Set the next event"                  );
43     pMenu->AddButton("         Quit        ",".q"         ,"close session"                       );
44   pMenu->Show();
45 }//Menu()
46 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
47 void General()
48 {         
49   TControlBar *pMenu = new TControlBar("vertical","General purpose",100,50);  
50     pMenu->AddButton("Debug ON","don();"                    ,"Switch debug on-off"                        );   
51     pMenu->AddButton("Debug OFF","doff();"                   ,"Switch debug on-off"                        );   
52     pMenu->AddButton("Geo GUI","geo();"                    ,"Shows geometry"                             ); 
53     pMenu->AddButton("Browser","new TBrowser;"             ,"Start ROOT TBrowser"                        );
54   pMenu->Show();  
55 }//General()
56 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
57 void SimData()
58 {
59   TControlBar *pMenu = new TControlBar("vertical","Sim data",310,50);  
60     pMenu->AddButton("Display ","hed();"    ,"Display Fast");
61     pMenu->AddButton("HITS QA"           ,"hqa()"     ,"QA plots for hits: hqa()");
62     pMenu->AddButton("Print stack"       ,"stack();"  ,"To print hits:     hp(evt)");
63     pMenu->AddButton("Print hits"        ,"hp();"     ,"To print hits:     hp(evt)");
64     pMenu->AddButton("Print sdigits"     ,"sp();"     ,"To print sdigits:  sp(evt)");
65     pMenu->AddButton("Print digits"      ,"dp();"     ,"To print digits:   dp(evt)");
66     pMenu->AddButton("Print clusters"    ,"cp();"     ,"To print clusters: cp(evt)");
67   pMenu->Show();         
68 }//SimData()
69 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
70 void RawData()
71 {
72   TControlBar *pMenu = new TControlBar("vertical","Raw data",580,50);  
73     pMenu->AddButton("ESD print"                       ,"ep();"                  ,"To print ESD info: ep()"         );  
74     pMenu->AddButton("ESD QA"                          ,"eq();"                  ,"To draw ESD hists: eq()"         );  
75     pMenu->AddButton("Clusters print"                  ,"cp();"                  ,"To print clusters: cp()"         );  
76     pMenu->AddButton("Clusters QA"                     ,"cq();"                  ,"To draw clusters hists: cq()"    );  
77     pMenu->AddButton("Print Matrix"                    ,"mp();"                  ,"To print prob matrix: mp()"      );  
78     pMenu->AddButton("Print occupancy"                 ,"r->OccupancyPrint(-1);" ,"To print occupancy"              );  
79     pMenu->AddButton("Print event summary  "           ,"r->SummaryOfEvent();"   ,"To print a summary of the event" );  
80   pMenu->Show();         
81 }//RawData()
82 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
83 void Test()
84 {         
85   TControlBar *pMenu = new TControlBar("vertical","Test",820,50);  
86     pMenu->AddButton("TEST Display "      ,"sed();"    ,"Display Fast");
87     pMenu->AddButton("Hits->Digits"       ,"thd();"                    ,"test hits->sdigits->digits"                 );   
88     pMenu->AddButton("Segmentation"       ,"ts()"                      ,"test segmentation methods"                  );
89     pMenu->AddButton("Test response"      ,"AliHMPIDParam::TestResp();","Test AliHMPIDParam response methods"         );
90     pMenu->AddButton("Print map"          ,"PrintMap();"               ,"Test AliHMPIDParam transformation methods"   );
91     pMenu->AddButton("Test Recon"         ,"rec();"                    ,"Test AliHMPIDRecon"                          );
92   pMenu->Show();  
93 }//Test()
94 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
95
96
97 void doff(){  Printf("DebugOFF");  AliLog::SetGlobalDebugLevel(0);}
98 void don() {  Printf("DebugON");   AliLog::SetGlobalDebugLevel(AliLog::kDebug);}
99
100 void geo (                       ) {gGeoManager->GetTopVolume()->Draw("ogl");}
101   
102 void du  (                       ) {h->Dump         (   );}                //utility display 
103
104 void hp  (                       ) {h->HitPrint  (nEvent);}   //print hits for requested event
105 void hq  (                       ) {h->HitQA     (   );}   //hits QA plots for all events 
106 void sp  (                       ) {h->SdiPrint  (nEvent);}   //print sdigits for requested event
107 void sq  (                       ) {h->SdiPrint  (nEvent);}   //print sdigits for requested event
108 void dp  (                       ) {h->DigPrint  (nEvent);}   //print digits for requested event
109 void dq  (                       ) {AliHMPIDReconstructor::DigQA     (al );}   //digits QA plots for all events
110 void cp  (                       ) {h->CluPrint  (nEvent);                 }   //print clusters for requested event
111 void cq  (                       ) {AliHMPIDReconstructor::CluQA     (al );}   //clusters QA plots for all events
112
113 void ep  (                       ) {AliHMPIDTracker::EsdQA(1);            } 
114 void eq  (                       ) {AliHMPIDTracker::EsdQA();             }                   
115 void mp  (Double_t probCut=0.7   ) {AliHMPIDTracker::MatrixPrint(probCut);}                   
116
117 void PrevEvent()                   {nEvent--;if(nEvent<0)nEvent=0;Printf("Current event is %i",nEvent);}
118 void NextEvent()                   {nEvent++;Printf("Current event is %i",nEvent);}
119 void stack(                     )  {AliHMPIDParam::Stack();}    
120 void tid  (Int_t tid,Int_t evt=0)  {AliHMPIDParam::Stack(evt,tid);} 
121 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
122 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
123 void PrintMap()
124 {
125  
126   Double_t r2d=TMath::RadToDeg();
127
128   Double_t x=AliHMPIDDigit::SizeAllX(),y=AliHMPIDDigit::SizeAllY();
129     
130   Printf("\n\n\n");                                       
131   
132   for(int ch=6;ch>=0;ch--){
133     AliHMPIDDigit dL(ch,0,1,0,0),dR(ch,0,1,67,0);
134     TVector3 lt=rp->Lors2Mars(ch,0,y);                                              TVector3 rt=rp->Lors2Mars(ch,x,y);
135                                        TVector3 ce=rp->Lors2Mars(ch,x/2,y/2);
136     TVector3 lb=rp->Lors2Mars(ch,0,0);                                              TVector3 rb=rp->Lors2Mars(ch,x,0);
137     
138     Printf(" ____________________________");                                       
139     Printf("|%6.2fcm            %6.2fcm|"         ,lt.Mag()                             , rt.Mag()       );
140     Printf("|%6.2fde            %6.2fde|"         ,lt.Theta()*r2d                       , rt.Theta()*r2d );
141     Printf("|%6.2fde            %6.2fde|"         ,lt.Phi()*r2d                         , rt.Phi()*r2d   );                                       
142     Printf("|                            |"                                                       );
143     Printf("|DDL %2i    %7.2fcm   DDL %2i|"       ,dL.DdlIdx()    ,  ce.Mag()           , dR.DdlIdx()    );
144     Printf("| 0x%x    %7.2fdeg   0x%x|"           ,dL.DdlId()     ,  ce.Theta()*r2d     , dR.DdlId()     );
145     Printf("|          %7.2fdeg        |"                         ,  ce.Phi()*r2d                        );
146     Printf("|                            |");                                                                              
147     Printf("|%6.2fcm            %6.2fcm|"         ,lb.Mag()                             , rb.Mag()       );
148     Printf("|%6.2fde            %6.2fde|"         ,lb.Theta()*r2d                       , rb.Theta()*r2d );
149     Printf("|%6.2fde     Ch%i    %6.2fde|"        ,lb.Phi()*r2d   ,  ch                 , rb.Phi()*r2d   );                                       
150     Printf(" ----------------------------");                                         
151   }
152   
153   Double_t m[3]; 
154   for(int i=0;i<1000;i++){
155     Float_t xout=0,xin=gRandom->Rndm()*130.60;
156     Float_t yout=0,yin=gRandom->Rndm()*126.16;
157     Int_t   c=gRandom->Rndm()*6;
158     rp->Lors2Mars(c,xin,yin,m);
159     rp->Mars2Lors(c,m,xout,yout);
160     if( (xin-xout) != 0) Printf("Problem in X");
161     if( (yin-yout) != 0) Printf("Problem in Y");
162   }                
163   
164   Int_t ddl,r,d,a,ch,raw,pc,px,py; AliHMPIDDigit dd;
165   
166   ddl=0;raw=0x2214000;r= 8;d=8;a=20;
167   ddl=1;raw=0x2214000;r= 8;d=8;a=20;
168   
169   
170   ddl=2;raw=0x08d6000;r= 2;d=3;a=22;
171   ddl=3;raw=0x08d6000;r= 2;d=3;a=22;
172   
173   
174   ddl=6;raw=0x592e000;r=22;d=4;a=46;ch=3;pc=4;px=55;py=5;dd.Raw(ddl,raw); 
175   Printf("(ch=%i,pc=%i,x=%2i,y=%2i) ddl=%i raw=0x%h (r=%2i,d=%2i,a=%2i)",
176            ch,   pc,  px,   py,     ddl,   raw,      r,    d,    a); dd.Print(); 
177   ddl=7;raw=0x592e000;r=22;d=4;a=46;ch=3;pc=1;
178   
179   
180 }//PrintMap()
181 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
182 void HitQA(Double_t cut,Double_t cutele,Double_t cutR)
183 {
184 // Provides a set of control plots intended primarily for charged particle flux analisys
185 // Arguments: cut (GeV)    - cut on momentum of any charged particles but electrons, 
186 //            cetele (GeV) - the same for electrons-positrons
187 //            cutR (cm)    - cut on production vertex radius (cylindrical system)        
188   gBenchmark->Start("HitsAna");
189   
190   Double_t cutPantiproton    =cut;
191   Double_t cutPkaonminus     =cut;
192   Double_t cutPpionminus     =cut;
193   Double_t cutPmuonminus     =cut;
194   Double_t cutPpositron      =cutele;
195                     
196   Double_t cutPelectron      =cutele;
197   Double_t cutPmuonplus      =cut;
198   Double_t cutPpionplus      =cut;
199   Double_t cutPkaonplus      =cut;
200   Double_t cutPproton        =cut;
201                        
202
203   TH2F *pEleHitRZ    =new TH2F("EleHitRZ"    ,Form("e^{+} e^{-} hit %s;z[cm];R[cm]" ,GetName())     , 400,-300,300 ,400,-500,500);   //R-z
204   TH2F *pEleHitRP    =new TH2F("EleHitRP"    ,Form("e^{+} e^{-} hit %s;p[GeV];R[cm]",GetName())     ,1000,-1  ,1   ,400,   0,550);   //R-p
205   TH1F *pEleAllP     =new TH1F("EleAllP"     ,     "e^{+} e^{-} all;p[GeV]"                         ,1000,-1  ,1                );  
206   TH1F *pEleHitP     =new TH1F("EleHitP"     ,Form("e^{+} e^{-} hit %s;p[GeV]"      ,GetName())     ,1000,-1  ,1                );   
207   TH1F *pMuoHitP     =new TH1F("MuoHitP"     ,Form("#mu^{-} #mu^{+} hit %s;p[GeV]"  ,GetName())     ,1000,-4  ,4                ); 
208   TH1F *pPioHitP     =new TH1F("PioHitP"     ,Form("#pi^{-} #pi^{+} hit %s;p[GeV]"  ,GetName())     ,1000,-4  ,4                ); 
209   TH1F *pKaoHitP     =new TH1F("KaoHitP"     ,Form("K^{-} K^{+} hit %s;p[GeV]"      ,GetName())     ,1000,-4  ,4                ); 
210   TH1F *pProHitP     =new TH1F("ProHitP"     ,Form("p^{-} p^{+} hit %s;p[GeV]"      ,GetName())     ,1000,-4  ,4                ); 
211   TH2F *pFlux        =new TH2F("flux"        ,Form("%s flux with Rvertex<%.1fcm"    ,GetName(),cutR),10  ,-5  ,5   , 10,0    ,10); //special text hist
212   TH2F *pVertex      =new TH2F("vertex"      ,Form("%s 2D vertex of HMPID hit;x;y"   ,GetName())     ,120 ,0   ,600 ,120,0    ,600); //special text hist
213   TH1F *pRho         =new TH1F("rho"         ,Form("%s r of HMPID hit"               ,GetName())     ,600 ,0   ,600); //special text hist
214   pFlux->SetStats(0);
215   pFlux->GetXaxis()->SetBinLabel(1 ,Form("p^{-}>%.3fGeV/c"   ,cutPantiproton));        
216   pFlux->GetXaxis()->SetBinLabel(2 ,Form("K^{-}>%.3fGeV/c"   ,cutPkaonminus ));        
217   pFlux->GetXaxis()->SetBinLabel(3 ,Form("#pi^{-}>%.3fGeV/c" ,cutPpionminus ));      
218   pFlux->GetXaxis()->SetBinLabel(4 ,Form("#mu^{-}>%.3fGeV/c" ,cutPmuonminus ));      
219   pFlux->GetXaxis()->SetBinLabel(5 ,Form("e^{+}>%.3fGeV/c"   ,cutPpositron  ));        
220   
221   pFlux->GetXaxis()->SetBinLabel(6 ,Form("e^{-}>%.3fGeV/c"   ,cutPelectron  ));        
222   pFlux->GetXaxis()->SetBinLabel(7 ,Form("#mu^{+}>%.3fGeV/c" ,cutPmuonplus  ));      
223   pFlux->GetXaxis()->SetBinLabel(8 ,Form("#pi^{+}>%.3fGeV/c" ,cutPpionplus  ));      
224   pFlux->GetXaxis()->SetBinLabel(9 ,Form("K^{+}>%.3fGeV/c"   ,cutPkaonplus  ));        
225   pFlux->GetXaxis()->SetBinLabel(10,Form("p^{+}>%.3fGeV/c"   ,cutPproton    ));        
226   
227   pFlux->GetYaxis()->SetBinLabel(1,"sum");  
228   pFlux->GetYaxis()->SetBinLabel(2,"ch1");  
229   pFlux->GetYaxis()->SetBinLabel(3,"ch2");  
230   pFlux->GetYaxis()->SetBinLabel(4,"ch3");  
231   pFlux->GetYaxis()->SetBinLabel(5,"ch4");  
232   pFlux->GetYaxis()->SetBinLabel(6,"ch5");  
233   pFlux->GetYaxis()->SetBinLabel(7,"ch6");  
234   pFlux->GetYaxis()->SetBinLabel(8,"ch7");  
235   pFlux->GetYaxis()->SetBinLabel(9,"prim"); 
236   pFlux->GetYaxis()->SetBinLabel(10,"tot");  
237   
238 //end of hists definition
239    
240   Int_t iNevents=fLoader->GetRunLoader()->GetAliRun()->GetEventsPerRun(),iCntPrimParts=0,iCntTotParts=0;
241 //load all needed trees   
242   fLoader->LoadHits(); 
243   fLoader->GetRunLoader()->LoadHeader(); 
244   fLoader->GetRunLoader()->LoadKinematics();  
245   
246   for(Int_t iEvtN=0;iEvtN < iNevents;iEvtN++){//events loop
247     fLoader->GetRunLoader()->GetEvent(iEvtN);
248     AliInfo(Form(" %i event processes",fLoader->GetRunLoader()->GetEventNumber()));
249     AliStack *pStack= fLoader->GetRunLoader()->Stack(); 
250     
251     for(Int_t iParticleN=0;iParticleN<pStack->GetNtrack();iParticleN++){//stack loop
252       TParticle *pPart=pStack->Particle(iParticleN);
253
254       if(iParticleN%10000==0) AliInfo(Form(" %i particles read",iParticleN));
255     
256       switch(pPart->GetPdgCode()){
257         case kProtonBar: pFlux->Fill(-4.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-4.5,8); break;
258         case kKMinus:    pFlux->Fill(-3.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-3.5,8); break;
259         case kPiMinus:   pFlux->Fill(-2.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-2.5,8); break;
260         case kMuonMinus: pFlux->Fill(-1.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-1.5,8); break;
261         case kPositron:  pFlux->Fill(-0.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-0.5,8); pEleAllP->Fill(-pPart->P()); break;
262       
263         case kElectron:  pFlux->Fill( 0.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 0.5,8); pEleAllP->Fill( pPart->P()); break;      
264         case kMuonPlus:  pFlux->Fill( 1.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 1.5,8); break;      
265         case kPiPlus:    pFlux->Fill( 2.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 2.5,8); break;      
266         case kKPlus:     pFlux->Fill( 3.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 3.5,8); break;      
267         case kProton:    pFlux->Fill( 4.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 4.5,8); break;            
268       }//switch
269     }//stack loop
270 //now hits analiser        
271     for(Int_t iEntryN=0;iEntryN < fLoader->TreeH()->GetEntries();iEntryN++){//TreeH loop
272       fLoader->TreeH()->GetEntry(iEntryN);                                  //get current entry (prim)                
273       for(Int_t iHitN=0;iHitN < Hits()->GetEntries();iHitN++){//hits loop
274         AliHMPIDHit *pHit = (AliHMPIDHit*)Hits()->At(iHitN);            //get current hit
275         TParticle  *pPart=pStack->Particle(pHit->GetTrack());      //get stack particle which produced the current hit
276         
277         if(pPart->GetPDG()->Charge()!=0&&pPart->Rho()>0.1) pVertex->Fill(pPart->Vx(),pPart->Vy()); //safe margin for sec.
278         if(pPart->GetPDG()->Charge()!=0) pRho->Fill(pPart->Rho()); //safe margin for sec.
279         if(pPart->R()>cutR) continue;                                   //cut on production radius (cylindrical system) 
280       
281         switch(pPart->GetPdgCode()){
282           case kProtonBar: if(pPart->P()>cutPantiproton) {pProHitP->Fill(-pPart->P()); pFlux->Fill(-4.5,pHit->Ch());}break;
283           case kKMinus   : if(pPart->P()>cutPkaonminus)  {pKaoHitP->Fill(-pPart->P()); pFlux->Fill(-3.5,pHit->Ch());}break;
284           case kPiMinus  : if(pPart->P()>cutPpionminus)  {pPioHitP->Fill(-pPart->P()); pFlux->Fill(-2.5,pHit->Ch());}break;
285           case kMuonMinus: if(pPart->P()>cutPmuonminus)  {pMuoHitP->Fill(-pPart->P()); pFlux->Fill(-1.5,pHit->Ch());}break;        
286           case kPositron : if(pPart->P()>cutPpositron)   {pEleHitP->Fill(-pPart->P()); pFlux->Fill(-0.5,pHit->Ch()); 
287                pEleHitRP->Fill(-pPart->P(),pPart->R());  pEleHitRZ->Fill(pPart->Vz(),pPart->R()); }break;
288           
289           case kElectron : if(pPart->P()>cutPelectron)   {pEleHitP->Fill( pPart->P()); pFlux->Fill( 0.5,pHit->Ch()); 
290                pEleHitRP->Fill( pPart->P(),pPart->R());  pEleHitRZ->Fill(pPart->Vz(),pPart->R()); }break;
291           case kMuonPlus : if(pPart->P()>cutPmuonplus)   {pMuoHitP->Fill( pPart->P()); pFlux->Fill( 1.5,pHit->Ch());}break;                     
292           case kPiPlus   : if(pPart->P()>cutPpionplus)   {pPioHitP->Fill( pPart->P()); pFlux->Fill( 2.5,pHit->Ch());}break;           
293           case kKPlus    : if(pPart->P()>cutPkaonplus)   {pKaoHitP->Fill( pPart->P()); pFlux->Fill( 3.5,pHit->Ch());}break;           
294           case kProton   : if(pPart->P()>cutPproton)     {pProHitP->Fill( pPart->P()); pFlux->Fill( 4.5,pHit->Ch());}break;
295         }
296       }//hits loop      
297     }//TreeH loop
298     iCntPrimParts +=pStack->GetNprimary();
299     iCntTotParts  +=pStack->GetNtrack();
300   }//events loop                        
301 //unload all loaded staff  
302   fLoader->UnloadHits();  
303   fLoader->GetRunLoader()->UnloadHeader(); 
304   fLoader->GetRunLoader()->UnloadKinematics();  
305 //Calculater some sums
306   Stat_t sum=0;
307 //sum row, sum over rows  
308   for(Int_t i=1;i<=pFlux->GetNbinsX();i++){
309     sum=0; for(Int_t j=2;j<=8;j++)    sum+=pFlux->GetBinContent(i,j);    
310     pFlux->SetBinContent(i,1,sum);
311   }
312     
313 //display everything  
314   new TCanvas("canvas1",Form("Events %i Nprims=%i Nparticles=%i",iNevents,iCntPrimParts,iCntTotParts),1000,900); pFlux->Draw("text");  gPad->SetGrid();  
315 //total prims and particles
316   TLatex latex; latex.SetTextSize(0.02);
317   sum=0; for(Int_t i=1;i<=pFlux->GetNbinsX();i++) sum+=pFlux->GetBinContent(i,10);    latex.DrawLatex(5.1,9.5,Form("%.0f",sum));
318   sum=0; for(Int_t i=1;i<=pFlux->GetNbinsX();i++) sum+=pFlux->GetBinContent(i,9);     latex.DrawLatex(5.1,8.5,Form("%.0f",sum));
319   for(Int_t iCh=0;iCh<7;iCh++) {
320     sum=0; for(Int_t i=1;i<=pFlux->GetNbinsX();i++) sum+=pFlux->GetBinContent(i,iCh+2);latex.DrawLatex(5.1,iCh+1.5,Form("%.0f",sum));
321   }  
322   sum=0; for(Int_t i=1;i<=pFlux->GetNbinsX();i++) sum+=pFlux->GetBinContent(i,1);    latex.DrawLatex(5.1,0.5,Form("%.0f",sum));
323   
324   new TCanvas("cEleAllP"   ,"e" ,200,100); pEleAllP->Draw();
325   new TCanvas("cEleHitRP"  ,"e" ,200,100); pEleHitRP->Draw();
326   new TCanvas("cEleHitRZ"  ,"e" ,200,100); pEleHitRZ->Draw();
327   new TCanvas("cEleHitP"   ,"e" ,200,100); pEleHitP->Draw();
328   new TCanvas("cMuoHitP"   ,"mu",200,100); pMuoHitP->Draw();
329   new TCanvas("cPioHitP"   ,"pi",200,100); pPioHitP->Draw();
330   new TCanvas("cKaoHitP"   ,"K" ,200,100); pKaoHitP->Draw();
331   new TCanvas("cProHitP"   ,"p" ,200,100); pProHitP->Draw();
332   new TCanvas("cVertex"    ,"2d vertex" ,200,100); pVertex->Draw();
333   new TCanvas("cRho"    ,"Rho of sec" ,200,100); pRho->Draw();
334   
335   gBenchmark->Show("HitsPlots");
336 }//HitQA()
337 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
338 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
339 void CluQA(AliRunLoader *pAL)
340 {
341 // Quality assesment plots for clusters. 
342 // This methode takes list of digits and form list of clusters again in order to 
343 // calculate cluster shape and cluster particle mixture    
344   AliLoader *pRL=pAL->GetDetectorLoader("HMPID");  AliHMPID *pRich=(AliHMPID*)pAL->GetAliRun()->GetDetector("HMPID");//get pointers for HMPID and HMPID loader
345   Int_t iNevt=pAL->GetNumberOfEvents();  if(iNevt==0)             {AliInfoClass("No events");return;}   
346                                          if(pRL->LoadDigits())    {AliInfoClass("No digits file");return;}
347                                             pAL->LoadHeader();
348                                             pAL->LoadKinematics(); 
349 //                                            AliStack *pStack=pAL->Stack();
350   TH1::AddDirectory(kFALSE);
351   
352         
353   TH1F*    pQ=new TH1F("HmpAllQ"  ,"Charge All"           ,4000 ,0  ,4000);// Q hists
354   TH1F* pCerQ=new TH1F("HmpCerQ"  ,"Charge Ckov"          ,4000 ,0  ,4000);
355   TH1F* pMipQ=new TH1F("HmpMipQ"  ,"Charge MIP"           ,4000 ,0  ,4000);
356   
357   TH1F*    pS=new TH1F("HmpCluSize"    ,"Cluster size;size"         ,100  ,0  ,100 );// size hists
358   TH1F* pCerS=new TH1F("HmpCluCerSize" ,"Ckov size;size"            ,100  ,0  ,100 );
359   TH1F* pMipS=new TH1F("HmpCluMipSize" ,"MIP size;size"             ,100  ,0  ,100 );
360   
361   TH2F*    pM=new TH2F("HmpCluMap"     ,"Cluster map;x [cm];y [cm]" ,1000 ,0  ,AliHMPIDDigit::SizeAllX(),1000,0,AliHMPIDDigit::SizeAllY()); // maps
362   TH2F* pMipM=new TH2F("HmpCluMipMap"  ,"MIP map;x [cm];y [cm]"     ,1000 ,0  ,AliHMPIDDigit::SizeAllX(),1000,0,AliHMPIDDigit::SizeAllY());
363   TH2F* pCerM=new TH2F("HmpCluCerMap"  ,"Ckov map;x [cm];y [cm]"    ,1000 ,0  ,AliHMPIDDigit::SizeAllX(),1000,0,AliHMPIDDigit::SizeAllY());
364  
365   
366   
367   for(Int_t iEvt=0;iEvt<iNevt;iEvt++){
368     pAL->GetEvent(iEvt);               
369     pRL->TreeD()->GetEntry(0); 
370     TClonesArray *pCluLst=new TClonesArray("AliHMPIDCluster");//tmp list of clusters for this event
371     
372     for(Int_t iCh=0;iCh<7;iCh++) AliHMPIDReconstructor::Dig2Clu(pRich->DigLst(iCh),pCluLst,kFALSE);//cluster finder for all chamber if any digits present
373     
374     for(Int_t iClu=0;iClu<pCluLst->GetEntriesFast();iClu++){
375       AliHMPIDCluster *pClu = (AliHMPIDCluster*)pCluLst->At(iClu);
376       Int_t cfm=0; for(Int_t iDig=0;iDig<pClu->Size();iDig++)  cfm+=pClu->Dig(iDig)->Ch(); //collect ckov-fee-mip structure of current cluster ?????
377       Int_t iNckov=cfm/1000000;      Int_t iNfee =cfm%1000000/1000;      Int_t iNmip =cfm%1000000%1000; 
378
379                                              pQ   ->Fill(pClu->Q()) ; pS   ->Fill(pClu->Size()) ; pM    ->Fill(pClu->X(),pClu->Y()); //all clusters  
380       if(iNckov!=0 && iNfee==0 && iNmip==0) {pCerQ->Fill(pClu->Q()) ; pCerS->Fill(pClu->Size()) ; pCerM ->Fill(pClu->X(),pClu->Y());}//ckov only cluster
381       if(iNckov==0 && iNfee==0 && iNmip!=0) {pMipQ->Fill(pClu->Q()) ; pMipS->Fill(pClu->Size()) ; pMipM ->Fill(pClu->X(),pClu->Y());}//mip only cluster
382                                        
383     }//clusters loop   
384     pCluLst->Clear();delete pCluLst;
385   }//events loop
386   
387   pRL->UnloadDigits(); pAL->UnloadKinematics(); pAL->UnloadHeader();
388   TCanvas *pC=new TCanvas("RichCluQA",Form("QA for cluster from %i events",iNevt),1000,900); pC->Divide(3,3);
389   pC->cd(1);    pM->Draw();          pC->cd(2);    pQ->Draw();       pC->cd(3);    pS->Draw();        
390   pC->cd(4); pMipM->Draw();          pC->cd(5); pMipQ->Draw();       pC->cd(6); pMipS->Draw();        
391   pC->cd(7); pCerM->Draw();          pC->cd(8); pCerQ->Draw();       pC->cd(9); pCerS->Draw();        
392 }//CluQA()
393 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
394 void OccupancyPrint(Int_t iEvtNreq)
395 {
396 //prints occupancy for each chamber in a given event
397   Int_t iEvtNmin,iEvtNmax;
398   if(iEvtNreq==-1){
399     iEvtNmin=0;
400     iEvtNmax=gAlice->GetEventsPerRun();
401   } else { 
402     iEvtNmin=iEvtNreq;iEvtNmax=iEvtNreq+1;
403   }
404     
405   if(GetLoader()->GetRunLoader()->LoadHeader()) return;    
406   if(GetLoader()->GetRunLoader()->LoadKinematics()) return;    
407   
408 //  Info("Occupancy","for event %i",iEvtN);
409   if(GetLoader()->LoadHits()) return;
410   if(GetLoader()->LoadDigits()) return;
411
412   
413   for(Int_t iEvtN=iEvtNmin;iEvtN<iEvtNmax;iEvtN++){    
414     Int_t nDigCh[7]={0,0,0,0,0,0,0};  
415     Int_t iChHits[7]={0,0,0,0,0,0,0};
416     Int_t nPrim[7]={0,0,0,0,0,0,0};
417     Int_t nSec[7]={0,0,0,0,0,0,0};
418     AliInfo(Form("events processed %i",iEvtN));
419     if(GetLoader()->GetRunLoader()->GetEvent(iEvtN)) return;    
420     AliStack *pStack = GetLoader()->GetRunLoader()->Stack();
421     for(Int_t iPrimN=0;iPrimN<GetLoader()->TreeH()->GetEntries();iPrimN++){//prims loop
422       GetLoader()->TreeH()->GetEntry(iPrimN);      
423       for(Int_t iHitN=0;iHitN<Hits()->GetEntries();iHitN++){
424         AliHMPIDHit *pHit = (AliHMPIDHit*)Hits()->At(iHitN);
425         if(pHit->E()>0){
426           iChHits[pHit->Ch()]++;
427           if(pStack->Particle(pHit->GetTrack())->Rho()<0.01) nPrim[pHit->Ch()]++;else nSec[pHit->Ch()]++;
428         }
429       }
430     }
431     
432     GetLoader()->TreeD()->GetEntry(0);
433     for(Int_t iCh=0;iCh<7;iCh++){
434       for(Int_t iDig=0;iDig<DigLst(iCh)->GetEntries();iDig++){
435         AliHMPIDDigit *pDig=(AliHMPIDDigit*)DigLst(iCh)->At(iDig);
436         nDigCh[pDig->Ch()]++;
437       }
438       Printf("Occupancy for chamber %i = %4.2f %% and charged prim tracks %i and sec. tracks %i with total %i",
439         iCh,Float_t(nDigCh[iCh])*100/AliHMPIDDigit::kPadAll,nPrim[iCh],nSec[iCh],iChHits[iCh]);
440     }      
441     
442     
443   }//events loop
444   GetLoader()->UnloadHits();
445   GetLoader()->UnloadDigits();
446   GetLoader()->GetRunLoader()->UnloadHeader();    
447   GetLoader()->GetRunLoader()->UnloadKinematics();    
448 }//OccupancyPrint()
449 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
450 void AliHMPID::SummaryOfEvent(Int_t iEvtN) const
451 {
452 //prints a summary for a given event
453   AliInfo(Form("Summary of event %i",iEvtN));
454   GetLoader()->GetRunLoader()->GetEvent(iEvtN);    
455   if(GetLoader()->GetRunLoader()->LoadHeader()) return;
456   if(GetLoader()->GetRunLoader()->LoadKinematics()) return;
457   AliStack *pStack=GetLoader()->GetRunLoader()->Stack();
458   
459   AliGenEventHeader* pGenHeader =  gAlice->GetHeader()->GenEventHeader();
460   if(pGenHeader->InheritsFrom("AliGenHijingEventHeader")) {
461     AliInfo(Form(" Hijing event with impact parameter b = %.2f (fm)",((AliGenHijingEventHeader*) pGenHeader)->ImpactParameter()));
462   }
463   Int_t nChargedPrimaries=0;
464   for(Int_t i=0;i<pStack->GetNtrack();i++) {
465     TParticle *pParticle = pStack->Particle(i);
466     if(pParticle->IsPrimary()&&pParticle->GetPDG()->Charge()!=0) nChargedPrimaries++;
467     }
468   AliInfo(Form("Total number of         primaries %i",pStack->GetNprimary()));
469   AliInfo(Form("Total number of charged primaries %i",nChargedPrimaries));
470   AliInfo(Form("Total n. of tracks in stack(+sec) %i",pStack->GetNtrack()));
471   GetLoader()->GetRunLoader()->UnloadHeader();
472   GetLoader()->GetRunLoader()->UnloadKinematics();
473 }//SummaryOfEvent()
474 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
475 void DrawZoom()
476 {
477 // Show info about current cursur position in status bar of the canvas
478 // Arguments: none
479 //   Returns: none     
480   TCanvas *pC=(TCanvas*)gPad; 
481   TRootCanvas *pRC= (TRootCanvas*)pC->GetCanvasImp();
482   TGStatusBar *pBar=pRC->GetStatusBar();
483   pBar->SetParts(5);
484   Float_t x=gPad->AbsPixeltoX(gPad->GetEventX());
485   Float_t y=gPad->AbsPixeltoY(gPad->GetEventY());
486   AliHMPIDDigit dig;dig.Manual1(1,x,y); UInt_t w32=0; 
487   if(IsInDead(x,y))
488     pBar->SetText("Out of sensitive area",4);    
489   else{
490     Int_t ddl=dig.Raw(w32);
491     pBar->SetText(Form("(p%i,x%i,y%i) ddl=%i 0x%x (r%i,d%i,a%i) (%.2f,%.2f)",
492         dig.Pc(),dig.PadPcX(),dig.PadPcY(),
493         ddl,w32,
494         dig.Row(),dig.Dilogic(),dig.Addr(),
495         dig.LorsX(),dig.LorsY()                            ),4);
496   }    
497   if(gPad->GetEvent()==1){
498     new TCanvas("zoom",Form("Row %i DILOGIC %i",dig.Row(),dig.Dilogic()));  
499   }
500 }//DrawZoom()
501 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
502 void DrawPc(Bool_t isFill) 
503
504 // Utility methode draws HMPID chamber PCs on event display.
505 // Arguments: none
506 //   Returns: none      
507 //  y6  ----------  ----------
508 //      |        |  |        |
509 //      |    4   |  |    5   |
510 //  y5  ----------  ----------
511 //
512 //  y4  ----------  ----------
513 //      |        |  |        |
514 //      |    2   |  |    3   |   view from electronics side
515 //  y3  ----------  ----------
516 //          
517 //  y2  ----------  ----------
518 //      |        |  |        |
519 //      |    0   |  |    1   |
520 //  y1  ----------  ----------
521 //      x1      x2  x3       x4
522   gPad->Range(-10,-10,AliHMPIDDigit::SizeAllX()+5,AliHMPIDDigit::SizeAllY()+5); 
523   Float_t x1=0,
524           x2=AliHMPIDDigit::SizePcX(),
525           x3=AliHMPIDDigit::SizePcX()+AliHMPIDDigit::SizeDead(),   
526           x4=AliHMPIDDigit::SizeAllX();
527   Float_t y1=0,
528           y2=  AliHMPIDDigit::SizePcY(),
529           y3=  AliHMPIDDigit::SizePcY()+AliHMPIDDigit::SizeDead(),
530           y4=2*AliHMPIDDigit::SizePcY()+AliHMPIDDigit::SizeDead(),
531           y5=  AliHMPIDDigit::SizeAllY()-AliHMPIDDigit::SizePcY(),
532           y6=  AliHMPIDDigit::SizeAllY();
533
534   Float_t xL[5]={x1,x1,x2,x2,x1}; //clockwise
535   Float_t xR[5]={x3,x3,x4,x4,x3};  
536   Float_t yD[5]={y1,y2,y2,y1,y1};
537   Float_t yC[5]={y3,y4,y4,y3,y3};  
538   Float_t yU[5]={y5,y6,y6,y5,y5};
539     
540   Float_t dX2=0.5*AliHMPIDDigit::SizePadX(),
541           dY2=0.5*AliHMPIDDigit::SizePadY() ;
542   
543   TLatex txt; txt.SetTextSize(0.01);
544   Int_t iColLeft=29,iColRight=41;
545   TPolyLine *pc=0;  TLine *pL;
546   AliHMPIDDigit dig;
547   for(Int_t iPc=0;iPc<AliHMPIDDigit::kPcAll;iPc++){
548     if(iPc==4) pc=new TPolyLine(5,xL,yU); if(iPc==5) pc=new TPolyLine(5,xR,yU); //draw PCs
549     if(iPc==2) pc=new TPolyLine(5,xL,yC); if(iPc==3) pc=new TPolyLine(5,xR,yC);
550     if(iPc==0) pc=new TPolyLine(5,xL,yD); if(iPc==1) pc=new TPolyLine(5,xR,yD);
551     (iPc%2)? pc->SetFillColor(iColLeft): pc->SetFillColor(iColRight);
552     if(!isFill){ pc->Draw();continue;}
553     
554     pc->Draw("f");
555     if(iPc%2) {dig.Manual2(0,iPc,79,25); txt.DrawText(dig.LorsX()+2,dig.LorsY(),Form("PC%i",dig.Pc()));}//print PC#    
556     
557     txt.SetTextAlign(32);
558     for(Int_t iRow=0;iRow<8 ;iRow++){//draw row lines (horizontal)
559       dig.Manual2(0,iPc,0,iRow*6);   //set digit to the left-down pad of this row
560       if(iPc%2) txt.DrawText(dig.LorsX()-1           ,dig.LorsY(),Form("%i",dig.PadPcY()));                  //print PadY#    
561                 txt.DrawText(dig.LorsX()-1+(iPc%2)*67,dig.LorsY()+2,Form("r%i",dig.Row()));                  //print Row#    
562       pL=new TLine(dig.LorsX()-dX2,dig.LorsY()-dY2,dig.LorsX()+AliHMPIDDigit::SizePcX()-dX2,dig.LorsY()-dY2);//draw horizontal line 
563       if(iRow!=0) pL->Draw(); 
564     }//row loop  
565     
566     txt.SetTextAlign(13);
567     for(Int_t iDil=0;iDil<10;iDil++){//draw dilogic lines (vertical)
568       dig.Manual2(0,iPc,iDil*8,0);       //set this digit to the left-down pad of this dilogic        
569                            txt.DrawText(dig.LorsX()  ,dig.LorsY()-1,Form("%i",dig.PadPcX()));                 //print PadX# 
570       if(iPc==4 || iPc==5) txt.DrawText(dig.LorsX()+2,dig.LorsY()+42,Form("d%i",dig.Dilogic()));              //print Dilogic#    
571       pL=new TLine(dig.LorsX()-dX2,dig.LorsY()-dY2,dig.LorsX()-dX2,dig.LorsY()+AliHMPIDDigit::SizePcY()-dY2); //draw vertical line
572       if(iDil!=0)pL->Draw();
573     }//dilogic loop        
574   }//PC loop      
575 }//DrawPc()
576 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
577 void hed()
578 {//event display from files
579   static TCanvas *pC=0;
580   static Int_t iEvt=0;
581   static Int_t iEvtTot=999;
582   static TFile *pEsdFl=0;
583   static TTree *pEsdTr=0;
584   static AliESD *pEsd=0;
585   if(!pC&&iEvt<iEvtTot){
586     if(hl==0) {Printf("hed: no HMPID loader");return;}
587     Printf("Opening session");
588     pEsdFl=TFile::Open("AliESDs.root");     if(!pEsdFl || !pEsdFl->IsOpen()) return;//open AliESDs.root
589     pEsdTr=(TTree*) pEsdFl->Get("esdTree"); if(!pEsdTr)                      return;//get ESD tree
590     pEsdTr->SetBranchAddress("ESD", &pEsd);
591     hl->LoadDigits(); hl->LoadRecPoints();
592     iEvtTot=pEsdTr->GetEntries();
593     pC=new TCanvas("hed","View from electronics side, IP is behind the picture.",1000,900);  pC->ToggleEventStatus(); pC->Divide(3,3);
594     pC->cd(7); TButton *pBtn=new TButton("Next","hed()",0,0,0.2,0.1);   pBtn->Draw(); 
595   }
596  
597   if(iEvt<iEvtTot){
598     pEsdTr->GetEntry(iEvt); al->GetEvent(iEvt); hl->TreeD()->GetEntry(0); hl->TreeR()->GetEntry(0);
599     TLatex txt; pC->cd(3);  gPad->Clear(); 
600     txt.SetTextSize(0.1);txt.DrawLatex(0.2,0.2,Form("Event %i (total %i)",iEvt,iEvtTot));
601     DrawEvt(pC,h->DigLst(),h->CluLst(),pEsd);
602     iEvt++;
603   }else{
604     Printf("--- No more events available...Bye.");
605     pC->Close();
606   }
607 }
608 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
609 void sed()
610 {
611
612   static TCanvas *pC1=0;
613   
614   if(!pC1){
615     pC1=new TCanvas("hed","Simulated evets-View from electronics side, IP is behind the picture.",1000,900); pC1->Divide(3,3);
616     pC1->cd(7); TButton *pBtn=new TButton("Next","sed()",0,0,0.2,0.1);   pBtn->Draw(); 
617   }
618
619
620   
621   AliHMPIDRecon rec;  
622
623   TClonesArray lh("AliHMPIDHit"); 
624   TClonesArray ls("AliHMPIDDigit"); 
625   TObjArray    ld(7); for(Int_t i=0;i<7;i++) ld.AddAt(new TClonesArray("AliHMPIDDigit"),i);
626   TObjArray    lc(7); for(Int_t i=0;i<7;i++) lc.AddAt(new TClonesArray("AliHMPIDCluster"),i);
627   AliESD esd;
628   
629   
630   EsdFromStack(&esd);
631   HitsFromEsd(&esd,&lh);
632   
633
634
635
636              AliHMPIDv1::Hit2Sdi(&lh,&ls);                               
637       AliHMPIDDigitizer::Sdi2Dig(&ls,&ld);     
638   AliHMPIDReconstructor::Dig2Clu(&ld,&lc);
639 //        AliHMPIDTracker::Recon(&esd,&cl);
640   
641   DrawEvt(pC1,&ld,&lc,&esd);  
642 }//SimEvt()
643 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
644 void DrawEvt(TCanvas *pC,TObjArray *pDigLst,TObjArray *pCluLst,AliESD *pEsd)
645 {//draws all the objects of current event
646
647   AliHMPIDRecon rec;  
648   TPolyMarker *pTxC[7];  TPolyMarker *pRin[7]; //intesections and rings
649   for(Int_t ch=0;ch<7;ch++){
650     pTxC[ch]=new TPolyMarker; pTxC[ch]->SetMarkerStyle(2); pTxC[ch]->SetMarkerColor(kRed); pTxC[ch]->SetMarkerSize(3);
651     pRin[ch]=new TPolyMarker; pRin[ch]->SetMarkerStyle(6); pRin[ch]->SetMarkerColor(kMagenta);
652   }
653   
654   
655   for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){//tracks loop to collect cerenkov rings and intersection points
656     AliESDtrack *pTrk=pEsd->GetTrack(iTrk);
657     Int_t ch=pTrk->GetHMPIDcluIdx();
658     if(ch<0) continue; //this track does not hit HMPID
659     ch/=1000000; 
660     Float_t th,ph,xPc,yPc; pTrk->GetHMPIDtrk(xPc,yPc,th,ph);  //get info on current track
661     pTxC[ch]->SetNextPoint(xPc,yPc);                          //add new intersection point
662     
663     Float_t ckov=pTrk->GetHMPIDsignal();  Float_t err=TMath::Sqrt(pTrk->GetHMPIDchi2());
664     
665     if(ckov>0){
666       rec.SetTrack(xPc,yPc,th,ph);
667      TVector2 pos;  for(int j=0;j<100;j++){rec.TracePhot(ckov,j*0.0628,pos); pRin[ch]->SetNextPoint(pos.X(),pos.Y());}      
668     }
669   }//tracks loop
670       
671   for(Int_t iCh=0;iCh<7;iCh++){//chambers loop    
672     switch(iCh){
673       case 6: pC->cd(1); break; case 5: pC->cd(2); break;
674       case 4: pC->cd(4); break; case 3: pC->cd(5); break; case 2: pC->cd(6); break;
675                                 case 1: pC->cd(8); break; case 0: pC->cd(9); break;
676     }
677     gPad->SetEditable(kTRUE); gPad->Clear();
678     DrawPc(0);
679     ((TClonesArray*)pDigLst->At(iCh))->Draw();  //draw digits
680     ((TClonesArray*)pCluLst->At(iCh))->Draw();  //draw clusters
681                             pTxC[iCh]->Draw();  //draw intersections
682                             pRin[iCh]->Draw();  //draw rings
683     gPad->SetEditable(kFALSE);
684   }//chambers loop
685 //  TLatex txt; txt.SetTextSize(0.02);
686 //  txt.DrawLatex(20,-5,Form("#theta=%.4f#pm%.5f with %2i #check{C}"          ,simCkov,simErr,simN));
687 //  txt.DrawLatex(25,-5,Form("#theta=%.4f#pm%.5f with %2i #check{C}"          ,recCkov,recErr,recN));
688 //  txt.DrawLatex(0 ,127,Form("#theta=%.2f#circ   #phi=%.2f#circ @(%.2f,%.2f) ",th*TMath::RadToDeg(),ph*TMath::RadToDeg(),radx,rady));
689 //  Printf("DIG------DIG---------DIG--------DIG------DIG------DIG");pDigLst->Print();Printf("");                   
690 }//DrawEvt()
691 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
692 void t1(Int_t case=1)
693 {
694   AliHMPIDDigit *d[10]; for(Int_t i=0;i<10;i++) d[i]=new AliHMPIDDigit;
695   
696   
697   Int_t iNdig;
698   
699   if(case==1){
700     iNdig=9;  
701   
702                                                               d[0]->Manual2(1,2,67,26, 33); 
703                                 d[1]->Manual2(1,2,66,25,431); d[2]->Manual2(1,2,67,25, 21);
704   d[3]->Manual2(1,2,65,24,127); d[4]->Manual2(1,2,66,24, 54); d[5]->Manual2(1,2,67,24,  5);
705   d[6]->Manual2(1,2,65,23, 20); d[7]->Manual2(1,2,66,23,  5); d[8]->Manual2(1,2,67,23,  6);
706   }else if(case==2){
707     iNdig=3;
708     d[0]->Manual2(0,0,36,14,  8); 
709     d[1]->Manual2(0,0,36,13, 33); d[2]->Manual2(0,0,37,13, 22);
710   }
711   
712   AliHMPIDCluster c;
713   for(int i=0;i<iNdig;i++) c.DigAdd(d[i]);  c.Print();
714   
715   
716   TClonesArray *cl=new TClonesArray("AliHMPIDCluster");
717   
718   c.Solve(cl,kTRUE);
719   Printf("");
720   
721   cl->Print();  
722 }//t1()
723 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
724 void EsdFromStack(AliESD *pEsd)
725 {
726   al->LoadHeader();al->LoadKinematics();
727   AliStack *pStk=al->Stack();  
728   
729   for(Int_t iTrk=0;iTrk<pStk->GetNtrack();iTrk++){//stack loop
730     TParticle *pPart=pStk->Particle(iTrk);
731     if(pPart->GetPDG()->Charge()==0) continue; //neutral particles are not reconstructed
732     if(pPart->GetFirstMother()>0)    continue; //do not consider secondaries
733     AliESDtrack trk(pPart);
734     pEsd->AddTrack(&trk);
735   }//stack loop  
736 }//EsdFromStack()
737 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
738 void HitsFromEsd(AliESD *pEsd, TClonesArray *pHitLst)
739 {
740   AliHMPIDRecon rec;
741   const Int_t kCerenkov=50000050,kFeedback=50000051;
742   Int_t hc=0; TVector2 pos;
743   for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){//tracks loop
744     AliESDtrack *pTrk=pEsd->GetTrack(iTrk);
745     Float_t xRa,yRa;
746     Int_t ch=AliHMPIDTracker::IntTrkCha(pTrk,xRa,yRa);
747     if(ch<0) continue; //this track does not hit HMPID
748     Float_t ckov=0.63;
749
750     Float_t th,ph,xPc,yPc,; pTrk->GetHMPIDtrk(xPc,yPc,th,ph);  
751                           new((*pHitLst)[hc++]) AliHMPIDHit(ch,200e-9,kProton  ,iTrk,xPc                ,yPc);                 //mip hit
752     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
753     for(int i=0;i<16;i++){
754       rec.TracePhot(ckov,gRandom->Rndm()*TMath::TwoPi(),pos);
755                           new((*pHitLst)[hc++]) AliHMPIDHit(ch,7.5e-9,kCerenkov,iTrk,pos.X(),pos.Y());}                      //photon hits  
756   }//tracks loop    
757 }
758 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++