RICH/menu.C

   1
   2 Int_t countContrib[7][3];
   3
   4 void ControlPlots()
   5 {
   6   Int_t iChamber=1;
   7
   8   TFile *pFile = new TFile("$(HOME)/plots.root","RECREATE");
   9   TH1F *pCqH1=new TH1F("ClusQ",   "Cluster Charge all chambers;q [QDC]",r->P()->MaxQdc(),0,r->P()->MaxQdc());
  10   TH1F *pCsH1=new TH1F("ClusSize","Cluster size all chambers;size [number of pads in cluster]",100,0,100);
  11   TH2F *pCmH2=new TH2F("ClusMap", "Cluster map;x [cm];y [cm]",1000,-r->P()->PcSizeX()/2,r->P()->PcSizeX()/2,
  12                                                              1000,-r->P()->PcSizeY()/2,r->P()->PcSizeY()/2);
  13
  14   TH1F *pCqMipH1=new TH1F("MipClusQ",   "MIP Cluster Charge all chambers;q [QDC]",r->P()->MaxQdc(),0,r->P()->MaxQdc());
  15   TH1F *pCsMipH1=new TH1F("MipClusSize","MIP Cluster size all chambers;size [number of pads in cluster]",100,0,100);
  16   TH2F *pCmMipH2=new TH2F("MipClusMap", "MIP Cluster map;x [cm];y [cm]",1000,-r->P()->PcSizeX()/2,r->P()->PcSizeX()/2,
  17                                                              1000,-r->P()->PcSizeY()/2,r->P()->PcSizeY()/2);
  18
  19   TH1F *pCqCerH1=new TH1F("CerClusQ",   "Cerenkov Cluster Charge all chambers;q [QDC]",r->P()->MaxQdc(),0,r->P()->MaxQdc());
  20   TH1F *pCsCerH1=new TH1F("CerClusSize","Cernekov Cluster size all chambers;size [number of pads in cluster]",100,0,100);
  21   TH2F *pCmCerH2=new TH2F("CerClusMap", "Cerenkov Cluster map;x [cm];y [cm]",1000,-r->P()->PcSizeX()/2,r->P()->PcSizeX()/2,
  22                                                              1000,-r->P()->PcSizeY()/2,r->P()->PcSizeY()/2);
  23   TH1F *pCqFeeH1=new TH1F("FeeClusQ",   "Feedback Cluster Charge all chambers;q [QDC]",r->P()->MaxQdc(),0,r->P()->MaxQdc());
  24   TH1F *pCsFeeH1=new TH1F("FeeClusSize","Feedback Cluster size all chambers;size [number of pads in cluster]",100,0,100);
  25   TH2F *pCmFeeH2=new TH2F("FeeClusMap", "Feedback Cluster map;x [cm];y [cm]",1000,-r->P()->PcSizeX()/2,r->P()->PcSizeX()/2,
  26                                                              1000,-r->P()->PcSizeY()/2,r->P()->PcSizeY()/2);
  27   Bool_t isClusters=!rl->LoadRecPoints();
  28   r->SetTreeAddress();
  29   for(Int_t iEventN=0;iEventN<a->GetEventsPerRun();iEventN++){//events loop
  30     al->GetEvent(iEventN);
  31     if(isClusters){
  32       rl->TreeR()->GetEntry(0);
  33       Int_t iTotalClusters=0;
  34       for(int i=1;i<=7;i++){//chambers loop
  35         iTotalClusters+=r->Clusters(i)->GetEntries();
  36         for(Int_t iClusterN=0;iClusterN<r->Clusters(i)->GetEntries();iClusterN++){//clusters loop
  37           AliRICHcluster *pClus=(AliRICHcluster*)r->Clusters(i)->At(iClusterN);
  38
  39           countContrib[i-1][0] += pClus->Nmips();
  40           countContrib[i-1][1] += pClus->Ncerenkovs();
  41           countContrib[i-1][2] += pClus->Nfeedbacks();
  42
  43           pCqH1->Fill(pClus->Q());
  44           pCsH1->Fill(pClus->Size());
  45           pCmH2->Fill(pClus->X(),pClus->Y());
  46
  47           if(pClus->IsPureMip()){ //Pure Mips
  48             pCqMipH1->Fill(pClus->Q());
  49             pCsMipH1->Fill(pClus->Size());
  50             pCmMipH2->Fill(pClus->X(),pClus->Y());
  51           }
  52           if(pClus->IsPureCerenkov()){ //Pure Photons
  53             pCqCerH1->Fill(pClus->Q());
  54             pCsCerH1->Fill(pClus->Size());
  55             pCmCerH2->Fill(pClus->X(),pClus->Y());
  56           }
  57           if(pClus->IsPureFeedback()){ //Pure Feedbacks
  58             pCqFeeH1->Fill(pClus->Q());
  59             pCsFeeH1->Fill(pClus->Size());
  60             pCmFeeH2->Fill(pClus->X(),pClus->Y());
  61           }
  62         }//clusters loop
  63       }//chambers loop
  64     }//isClusters
  65     Info("ControlPlots","Event %i processed.",iEventN);
  66   }//events loop
  67   if(isClusters) rl->UnloadRecPoints();
  68
  69   pFile->Write();
  70   delete pFile;
  71   for(Int_t i=0;i<7;i++)
  72     cout <<" chamber " << i+1 << " n. mips " << countContrib[i][0] << " n. ckovs " << countContrib[i][1] << " n. fdbks " << countContrib[i][2] << endl;
  73 }//void ControlPlots()
  74 //__________________________________________________________________________________________________
  75 void MainTrank()
  76 {
  77   TStopwatch sw;TDatime time;
  78   OLD_S_SD(); SD_D();   AliRICHClusterFinder *z=new AliRICHClusterFinder(r); z->Exec();//delete z;
  79   cout<<"\nInfo in <MainTrank>: Start time: ";time.Print();
  80     cout<<"Info in <MainTrank>: Stop  time: ";time.Set();  time.Print();
  81     cout<<"Info in <MainTrank>: Time  used: ";sw.Print();
  82 }
  83 //__________________________________________________________________________________________________
  84 void sh()
  85 {
  86   if(rl->LoadHits()) return;
  87
  88   Int_t iTotalHits=0;
  89   for(Int_t iPrimN=0;iPrimN<rl->TreeH()->GetEntries();iPrimN++){//prims loop
  90     rl->TreeH()->GetEntry(iPrimN);
  91     r->Hits()->Print();
  92     iTotalHits+=r->Hits()->GetEntries();
  93   }
  94   Info("sh","totally %i hits",iTotalHits);
  95   rl->UnloadHits();
  96 }
  97 //__________________________________________________________________________________________________
  98 void ssp()
  99 {
 100   if(rl->LoadHits()) return;
 101
 102   Int_t iTotalSpecials=0;
 103   for(Int_t iPrimN=0;iPrimN<rl->TreeH()->GetEntries();iPrimN++){//prims loop
 104     rl->TreeH()->GetEntry(iPrimN);
 105     r->Specials()->Print();
 106     iTotalSpecials+=r->Specials()->GetEntries();
 107   }
 108   Info("ssp","totally %i specials",iTotalSpecials);
 109   rl->UnloadHits();
 110 }
 111 //__________________________________________________________________________________________________
 112 void ss()
 113 {
 114   if(rl->LoadSDigits()) return;
 115   rl->TreeS()->GetEntry(0);
 116   r->SDigits()->Print();
 117   Info("ss","totally %i sdigits",r->SDigits()->GetEntries());
 118   rl->UnloadSDigits();
 119 }
 120 //__________________________________________________________________________________________________
 121 void sd()
 122 {
 123   if(rl->LoadDigits()) return;
 124   rl->TreeD()->GetEntry(0);
 125   Int_t iTotalDigits=0;
 126   for(int i=1;i<=7;i++){
 127     r->Digits(i)->Print();
 128     iTotalDigits+=r->Digits(i)->GetEntries();
 129   }
 130   Info("sd","totally %i digits",iTotalDigits);
 131   rl->UnloadDigits();
 132 }
 133
 134 void sc()
 135 {
 136   if(rl->LoadRecPoints()) return;
 137   r->SetTreeAddress();
 138   rl->TreeR()->GetEntry(0);
 139   for(int i=1;i<=7;i++) r->Clusters(i)->Print();
 140   rl->UnloadRecPoints();
 141 }
 142
 143 Double_t r2d = TMath::RadToDeg();
 144 Double_t d2r = TMath::DegToRad();
 145
 146 void DisplFast(){ AliRICHDisplFast *d = new AliRICHDisplFast();  d->Exec();}
 147
 148
 149 void C_R()
 150 {
 151   AliRICHRecon *detect = new AliRICHRecon("RICH patrec algorithm","Reconstruction");
 152
 153
 154   for (int nev=0; nev< a->GetEventsPerRun(); nev++) {    // Event Loop
 155     al->GetEvent(nev);
 156     cout <<endl<< "Processing event:" <<nev<<endl;
 157     detect->StartProcessEvent();
 158   } // event loop
 159   delete detect;
 160 }
 161 //__________________________________________________________________________________________________
 162 void D_C()
 163 {
 164   TStopwatch sw;TDatime time;
 165
 166    AliRICHClusterFinder *z=new AliRICHClusterFinder(r); z->Exec();
 167
 168    cout << endl;
 169    cout << "Info in Digits->Clusters: Start time: ";time.Print();
 170    cout << "Info in Digits->Clusters: Stop  time: ";time.Set();  time.Print();
 171    cout << "Info in Digits->Clusters: Time  used: ";sw.Print();
 172 }
 173 //__________________________________________________________________________________________________
 174 void OLD_S_SD()
 175 {
 176   Info("OLD_S_SD","Start.");
 177   rl->LoadHits();
 178   for(Int_t iEventN=0;iEventN<a->GetEventsPerRun();iEventN++){//events loop
 179     al->GetEvent(iEventN);    Info("OLD_S_SD","Processing event %i",iEventN);
 180
 181     rl->MakeTree("S");  r->MakeBranch("S");
 182     r->ResetSDigits();  r->ResetSpecialsOld();
 183
 184     for(Int_t iPrimN=0;iPrimN<rl->TreeH()->GetEntries();iPrimN++){//prims loop
 185       rl->TreeH()->GetEntry(iPrimN);
 186       for(Int_t i=0;i<r->Specials()->GetEntries();i++){//specials loop
 187         Int_t padx=1+ ((AliRICHSDigit*)r->Specials()->At(i))->PadX()+r->Param()->NpadsX()/2;
 188         Int_t pady=1+ ((AliRICHSDigit*)r->Specials()->At(i))->PadY()+r->Param()->NpadsY()/2;
 189         Double_t q=  ((AliRICHSDigit*)r->Specials()->At(i))->QPad();
 190
 191         Int_t hitN= ((AliRICHSDigit*)r->Specials()->At(i))->HitNumber()-1;//!!! important -1
 192         Int_t chamber=((AliRICHhit*)r->Hits()->At(hitN))->C();
 193         Int_t tid=((AliRICHhit*)r->Hits()->At(hitN))->GetTrack();
 194         Int_t pid=((AliRICHhit*)r->Hits()->At(hitN))->Pid();
 195         if(padx<1 || padx>r->Param()->NpadsX() ||pady<1 || pady>r->Param()->NpadsY())
 196           Warning("OLD_S_SD","pad is out of valid range padx= %i pady=%i event %i",padx,pady,iEventN);
 197         else
 198           r->AddSDigit(chamber,padx,pady,q,pid,tid);
 199       }//specials loop
 200     }//prims loop
 201     rl->TreeS()->Fill();
 202     rl->WriteSDigits("OVERWRITE");
 203   }//events loop
 204     rl->UnloadHits();     rl->UnloadSDigits();
 205   Info("OLD_S_SD","Stop.");
 206 }//OLD_S_SD()
 207 //__________________________________________________________________________________________________
 208 void H_SD()
 209 {
 210   Info("H_SD","Start.");
 211
 212   for(Int_t iEventN=0;iEventN<a->GetEventsPerRun();iEventN++){//events loop
 213     al->GetEvent(iEventN);
 214
 215     if(!rl->TreeH()) rl->LoadHits();//from
 216     if(!rl->TreeS()) rl->MakeTree("S");    r->MakeBranch("S");//to
 217
 218     for(Int_t iPrimN=0;iPrimN<rl->TreeH()->GetEntries();iPrimN++){//prims loop
 219       rl->TreeH()->GetEntry(iPrimN);
 220       for(Int_t iHitN=0;iHitN<3;iHitN++){//hits loop  ???
 221         AliRICHhit *pHit=r->Hits()->At(iHitN);
 222         TVector3 globX3(pHit->X(),pHit->Y(),pHit->Z());
 223         TVector3 locX3=r->C(pHit->C())->Glob2Loc(globX3);
 224
 225         Int_t sector;
 226         Int_t iTotQdc=r->Param()->Loc2TotQdc(locX3,pHit->Eloss(),pHit->Pid(),sector);
 227
 228         Int_t iPadXmin,iPadXmax,iPadYmin,iPadYmax;
 229         r->Param()->Loc2Area(locX3,iPadXmin,iPadYmin,iPadXmax,iPadYmax);
 230         cout<<"left-down=("<<iPadXmin<<","<<iPadYmin<<") right-up=("<<iPadXmax<<','<<iPadYmax<<')'<<endl;
 231         for(Int_t iPadY=iPadYmin;iPadY<=iPadYmax;iPadY++)
 232           for(Int_t iPadX=iPadXmin;iPadX<=iPadXmax;iPadX++){
 233             Double_t padQdc=iTotQdc*r->Param()->Loc2PadFrac(locX3,iPadX,iPadY);
 234             if(padQdc>0.1)r->AddSDigit(pHit->C(),iPadX,iPadY,padQdc,pHit->GetTrack());
 235           }
 236       }//hits loop
 237     }//prims loop
 238     rl->TreeS()->Fill();
 239     rl->WriteSDigits("OVERWRITE");
 240   }//events loop
 241
 242   rl->UnloadHits();
 243   rl->UnloadSDigits();
 244   Info("H_SD","Stop.");
 245 }//H_SD()
 246 //__________________________________________________________________________________________________
 247 void SD_D()
 248 {
 249   Info("SD_D","Start.");
 250   extern Int_t kBad;
 251   r->Param()->GenSigmaThMap();
 252   rl->LoadSDigits();
 253
 254   for(Int_t iEventN=0;iEventN<a->GetEventsPerRun();iEventN++){//events loop
 255     al->GetEvent(iEventN);    cout<<"Event "<<iEventN<<endl;
 256     rl->MakeTree("D");r->MakeBranch("D"); //create TreeD with RICH branches
 257     r->ResetSDigits();r->ResetDigits();//reset lists of sdigits and digits
 258     rl->TreeS()->GetEntry(0);
 259     r->SDigits()->Sort();
 260
 261     Int_t combiPid,chamber,x,y,tid[3],id; Double_t q;
 262     Int_t iNdigitsPerPad;//how many sdigits for a given pad
 263     const int kBad=-101;//??? to be removed in code
 264     for(Int_t i=0;i<r->SDigits()->GetEntries();i++){//sdigits loop (sorted)
 265       AliRICHdigit *pSdig=(AliRICHdigit*)r->SDigits()->At(i);
 266       if(pSdig->Id()==id){//still the same pad
 267         iNdigitsPerPad++;
 268         q+=pSdig->Q();
 269         combiPid+=pSdig->CombiPid();
 270         if(iNdigitsPerPad<=3)
 271           tid[iNdigitsPerPad-1]=pSdig->Tid(0);
 272 //        else
 273 //          Info("","More then 3 sdigits for the given pad");
 274       }else{//new pad, add the pevious one
 275         if(id!=kBad&&r->Param()->IsOverTh(chamber,x,y,q)) {
 276            r->AddDigit(chamber,x,y,q,combiPid,tid);
 277          }
 278         combiPid=pSdig->CombiPid();chamber=pSdig->C();id=pSdig->Id();
 279         x=pSdig->X();y=pSdig->Y();
 280         q=pSdig->Q();
 281         tid[0]=pSdig->Tid(0);
 282         iNdigitsPerPad=1;tid[1]=tid[2]=kBad;
 283       }
 284     }//sdigits loop (sorted)
 285
 286     if(r->SDigits()->GetEntries()&&r->Param()->IsOverTh(chamber,x,y,q))
 287       r->AddDigit(chamber,x,y,q,combiPid,tid);//add the last digit
 288
 289     rl->TreeD()->Fill();
 290     rl->WriteDigits("OVERWRITE");
 291   }//events loop
 292   rl->UnloadSDigits();     rl->UnloadDigits();
 293   r->ResetSDigits();r->ResetDigits();//reset lists of sdigits and digits
 294   Info("SD_D","Stop.");
 295 }//SD_D()
 296 //__________________________________________________________________________________________________
 297 void Show()
 298 {
 299   cout<<endl;
 300   al->LoadHeader();  al->LoadKinematics();
 301
 302   rl->LoadHits();
 303     Bool_t isSdigits=!rl->LoadSDigits();
 304       Bool_t isClusters=!rl->LoadRecPoints();
 305         Bool_t isDigits=!rl->LoadDigits();//loaders
 306
 307   cout<<endl;  cout<<endl;
 308   for(Int_t iEventN=0;iEventN<a->GetEventsPerRun();iEventN++){//events loop
 309     Int_t iNparticles=a->GetEvent(iEventN);
 310     Int_t iNprims=rl->TreeH()->GetEntries();
 311
 312     Int_t iTotalHits=0,iTotalCerenkovs=0,iTotalSpecials=0;
 313     for(Int_t iPrimN=0;iPrimN<iNprims;iPrimN++){//prims loop
 314       rl->TreeH()->GetEntry(iPrimN);
 315       iTotalHits+=r->Hits()->GetEntries();
 316       iTotalCerenkovs+=r->Cerenkovs()->GetEntries();
 317       iTotalSpecials+=r->Specials()->GetEntries();
 318       TParticle *pPrim=al->Stack()->Particle(iPrimN);
 319       Info("Show","Evt %4i prim %4i has %4i hits %5i cerenkovs and %5i specials from %s (,%7.2f,%7.2f)",
 320                            iEventN,
 321                                     iPrimN,
 322                                              r->Hits()->GetEntries(),
 323                                                       r->Cerenkovs()->GetEntries(),
 324                                                                         r->Specials()->GetEntries(),
 325                                                                                          pPrim->GetName(),
 326                                                                                  pPrim->Theta()*r2d,pPrim->Phi()*r2d);
 327     }//prims loop
 328     Info("Show-HITS","Evt %i total:  %i particles %i primaries %i hits %i cerenkovs %i specials",
 329                         iEventN,   iNparticles, iNprims,     iTotalHits,iTotalCerenkovs,iTotalSpecials);
 330     if(isSdigits){
 331       rl->TreeS()->GetEntry(0);
 332       Info("Show-SDIGITS","Evt %i contains %5i sdigits",iEventN,r->SDigits()->GetEntries());
 333     }
 334     if(isDigits){
 335       rl->TreeD()->GetEntry(0);
 336       for(int i=1;i<=7;i++)
 337         Info("Show-DIGITS","Evt %i chamber %i contains %5i digits",
 338                                  iEventN,   i,           r->Digits(i)->GetEntries());
 339     }
 340     if(isClusters){
 341       rl->TreeR()->GetEntry(0);
 342       for(int i=1;i<=7;i++)
 343         Info("Show-CLUSTERS","Evt %i chamber %i contains %5i clusters",
 344                                  iEventN,   i,           r->Clusters(i)->GetEntries());
 345     }
 346     cout<<endl;
 347   }//events loop
 348   rl->UnloadHits();
 349     if(isSdigits) rl->UnloadSDigits();
 350       if(isDigits) rl->UnloadDigits();
 351         if(isClusters) rl->UnloadRecPoints();
 352   al->UnloadHeader();
 353   al->UnloadKinematics();
 354   cout<<endl;
 355 }//void Show()
 356 //__________________________________________________________________________________________________
 357
 358
 359 AliRun *a;
 360 AliRunLoader *al;
 361 AliLoader *rl,*tl,*il;
 362
 363 AliRICH *r;
 364
 365 Bool_t ReadAlice()
 366 {
 367   Info("ReadAlice","Tring to read ALICE from SIMULATED FILE.");
 368   AliLoader::SetDebug(0);
 369   if(gAlice) delete gAlice;
 370   if(!(al=AliRunLoader::Open("galice.root","AlicE","update"))){
 371     gSystem->Exec("rm -rf *.root *.dat");
 372     Error("ReadAlice","galice.root broken, removing all this garbage then init new one");
 373     new AliRun("gAlice","Alice experiment system");
 374     gAlice->SetDebug(-1);
 375     gAlice->Init("ConfigRich.C");
 376     r=(AliRICH*)gAlice->GetDetector("RICH");
 377     return kFALSE;
 378   }
 379   al->LoadgAlice();
 380   if(!gAlice) Fatal("ReadAlice","No gAlice in file");
 381   a=al->GetAliRun();
 382   a->SetDebug(0);
 383 //RICH
 384   if(!(r=(AliRICH*)gAlice->GetDetector("RICH"))) Warning("RICH/menu.C::ReadAlice","No RICH in file");
 385   r->SetDebug(0);
 386   if(!(rl=al->GetLoader("RICHLoader")))          Warning("RICH/menu.C::ReadAlice","No RICH loader in file");
 387
 388   Info("ReadAlice","Run contains %i event(s)",gAlice->GetEventsPerRun());
 389   return kTRUE;
 390 }
 391 //__________________________________________________________________________________________________
 392 void PrintGeo(Float_t rotDeg=0)
 393 {
 394   AliRICHParam *p=new AliRICHParam;
 395   Double_t r=p->Offset();
 396   Double_t kP=p->AngleXY();
 397   Double_t kT=p->AngleYZ();
 398   Double_t kRot;
 399
 400   if(rotDeg==0)
 401     kRot=p->AngleRot();
 402   else
 403     kRot=rotDeg*deg;
 404
 405   cout<<endl;
 406   Double_t  phi=90*deg+kRot+kP,theta=90*deg+kT;
 407   Info("   menu for          1","r=%8.3f theta=%5.1f phi=%5.1f x=%8.3f y=%8.3f z=%8.3f",
 408                                  r,      theta*r2d,  phi*r2d,
 409                                                                r*sin(theta)*cos(phi),
 410                                                                        r*sin(theta)*sin(phi),
 411                                                                                r*cos(theta));
 412
 413   phi=90*deg+kRot+kP,theta=90*deg;
 414   Info("   menu for          2","r=%8.3f theta=%5.1f phi=%5.1f x=%8.3f y=%8.3f z=%8.3f",
 415                                  r,      theta*r2d,  phi*r2d,
 416                                                                r*sin(theta)*cos(phi),
 417                                                                        r*sin(theta)*sin(phi),
 418                                                                                r*cos(theta));
 419
 420   phi=90*deg+kRot,theta=90*deg-kT;
 421   Info("   menu for          3","r=%8.3f theta=%5.1f phi=%5.1f x=%8.3f y=%8.3f z=%8.3f",
 422                                  r,      theta*r2d,  phi*r2d,
 423                                                                r*sin(theta)*cos(phi),
 424                                                                        r*sin(theta)*sin(phi),
 425                                                                                r*cos(theta));
 426
 427
 428   phi=90*deg+kRot,theta=90*deg;
 429   Info("   menu for          4","r=%8.3f theta=%5.1f phi=%5.1f x=%8.3f y=%8.3f z=%8.3f",
 430                                  r,      theta*r2d,  phi*r2d,
 431                                                                r*sin(theta)*cos(phi),
 432                                                                        r*sin(theta)*sin(phi),
 433                                                                                r*cos(theta));
 434
 435
 436   phi=90*deg+kRot,theta=90*deg+kT;
 437   Info("   menu for          5","r=%8.3f theta=%5.1f phi=%5.1f x=%8.3f y=%8.3f z=%8.3f",
 438                                  r,      theta*r2d,  phi*r2d,
 439                                                                 r*sin(theta)*cos(phi),
 440                                                                        r*sin(theta)*sin(phi),
 441                                                                                r*cos(theta));
 442
 443
 444   phi=90*deg+kRot-kP,theta=90*deg;
 445   Info("   menu for          6","r=%8.3f theta=%5.1f phi=%5.1f x=%8.3f y=%8.3f z=%8.3f",
 446                                  r,      theta*r2d,  phi*r2d,
 447                                                                r*sin(theta)*cos(phi),
 448                                                                        r*sin(theta)*sin(phi),
 449                                                                                r*cos(theta));
 450
 451   phi=90*deg+kRot-kP,theta=90*deg+kT;
 452   Info("   menu for          7","r=%8.3f theta=%5.1f phi=%5.1f x=%8.3f y=%8.3f z=%8.3f",
 453                                  r,      theta*r2d,  phi*r2d,
 454                                                                r*sin(theta)*cos(phi),
 455                                                                        r*sin(theta)*sin(phi),
 456                                                                                r*cos(theta));
 457
 458   delete p;
 459 }//PrintGeo()
 460 //__________________________________________________________________________________________________
 461 void TestResponse()
 462 {
 463   TCanvas *pC=new TCanvas("c","Amplification test",900,800);
 464   pC->Divide(1,2);
 465   pC->cd(1);
 466   TF1 *pF1=new TF1("f1","9e-6*pow(x,4)+2e-7*pow(x,3)-0.0316*pow(x,2)-3e-4*x+25.367",-70,70);
 467   pF1->Draw();
 468
 469   pC->cd(2);
 470
 471   const Int_t nPoints=8;
 472   THStack *pStack=new THStack("stack","photons");
 473   TLegend *pLeg=new TLegend(0.6,0.2,0.9,0.5,"legend");
 474   TH1F *apH[nPoints];
 475
 476   Double_t starty=AliRICHParam::DeadZone()/2;
 477   Double_t deltay=AliRICHParam::SectorSizeY()/nPoints;
 478
 479   for(int i=0;i<nPoints;i++){
 480     apH[i]=new TH1F(Form("h%i",i),"Qdc for Photon;QDC;Counts",1000,0,1000); apH[i]->SetLineColor(i);
 481     pStack->Add(apH[i]);
 482     pLeg->AddEntry(apH[i],Form("@(0,%5.2f->%5.2f)",starty+i*deltay,starty+i*deltay-AliRICHParam::SectorSizeY()/2));
 483   }
 484
 485
 486   TVector3 x3(0,0,0);
 487   Int_t sector=10;
 488 //  AliRICHParam::ResetWireSag();
 489   for(Int_t i=0;i<10000;i++){//events loop
 490     for(int j=0;j<nPoints;j++){
 491       x3.SetY(starty-j*deltay);
 492       apH[j]->Fill(AliRICHParam::Loc2TotQdc(x3,400e-9,500000,sector));
 493     }
 494   }
 495   pStack->Draw("nostack");
 496   pLeg->Draw();
 497 }//TestResponse()
 498 //__________________________________________________________________________________________________
 499 void TestSD()
 500 {
 501   Info("TestSD","Creating test sdigits.");
 502   rl->MakeTree("S");r->MakeBranch("S");
 503
 504     for(Int_t iEventN=0;iEventN<a->GetEventsPerRun();iEventN++){//events loop
 505     al->GetEvent(iEventN);
 506
 507     if(!rl->TreeH()) rl->LoadHits();//from
 508     if(!rl->TreeS()) rl->MakeTree("S");    r->MakeBranch("S");//to
 509
 510     for(Int_t iPrimN=0;iPrimN<rl->TreeH()->GetEntries();iPrimN++){//prims loop
 511       rl->TreeH()->GetEntry(iPrimN);
 512       for(Int_t iHitN=0;iHitN<3;iHitN++){//hits loop  ???
 513         AliRICHhit *pHit=r->Hits()->At(iHitN);
 514         TVector3 globX3(pHit->X(),pHit->Y(),pHit->Z());
 515         TVector3 locX3=r->C(pHit->C())->Glob2Loc(globX3);
 516
 517         Int_t sector;
 518         Int_t iTotQdc=r->Param()->Loc2TotQdc(locX3,pHit->Eloss(),pHit->Pid(),sector);
 519
 520         Int_t iPadXmin,iPadXmax,iPadYmin,iPadYmax;
 521         r->Param()->Loc2Area(locX3,iPadXmin,iPadYmin,iPadXmax,iPadYmax);
 522         cout<<"left-down=("<<iPadXmin<<","<<iPadYmin<<") right-up=("<<iPadXmax<<','<<iPadYmax<<')'<<endl;
 523         for(Int_t iPadY=iPadYmin;iPadY<=iPadYmax;iPadY++)
 524           for(Int_t iPadX=iPadXmin;iPadX<=iPadXmax;iPadX++){
 525             Float_t iPadQdc=iTotQdc*r->Param()->Loc2PadFrac(locX3,iPadX,iPadY);
 526             Int_t padx,pady; r->Param()->Loc2Pad(locX3.X(),locX3.Y(),padx,pady);
 527             cout<<"hit="<<iHitN<<" ("<<locX3.X()<<','<<locX3.Y()<<")("<<padx<<','<<pady<<") cur pad("<<iPadX<<","<<iPadY<<") qtot="<<iTotQdc<<" qfrac="<<r->Param()->Loc2PadFrac(locX3,iPadX,iPadY)<<endl;
 528           }
 529 //            r->AddSDigit(pHit->C(),padx,pady,r->Param()->Local2PadQdc(localX3,padx,pady),pHit->GetTrack());
 530       }//hits loop
 531     }//prims loop
 532     rl->TreeS()->Fill();
 533     rl->WriteSDigits("OVERWRITE");
 534   }//events loop
 535
 536   rl->UnloadHits();
 537   rl->UnloadSDigits();
 538
 539   rl->TreeS()->Fill();
 540   rl->WriteSDigits("OVERWRITE");
 541   rl->UnloadSDigits();
 542   cout<<endl;r->Sdigits()->Print();
 543   r->ResetSDigits();
 544   Info("TestSdigits","Stop.");
 545 }//void TestSdigits()
 546 //__________________________________________________________________________________________________
 547 void TestC()
 548 {
 549   Info("TestC","Creating test clusters.");
 550   rl->MakeTree("R");r->MakeBranch("R");
 551
 552   AliRICHcluster c;
 553   c.AddDigit(new AliRICHdigit(1,20,21,200,1,2,3));
 554   c.AddDigit(new AliRICHdigit(1,22,21,250,1,2,3));
 555   c.CoG();
 556
 557   r->AddCluster(c);
 558
 559   rl->TreeR()->Fill();
 560   rl->WriteRecPoints("OVERWRITE");
 561   rl->UnloadRecPoints();
 562   r->ResetClusters();
 563
 564   Info("TestC","Stop.");
 565 }//TestC()
 566 //__________________________________________________________________________________________________
 567 void TestSeg()
 568 {
 569   AliRICHParam *p=r->Param();
 570   Int_t padx,pady;
 571   Double_t x,y;
 572   Double_t dz=p->DeadZone();
 573   Double_t sx=p->SectorSizeX(); Double_t sy=p->SectorSizeY();  Double_t px=p->PcSizeX(); Double_t py=p->PcSizeY();
 574   cout<<endl;
 575   Info("  1-  1","sec=%i padx=%3i pady=%3i",p->Loc2Pad(-px/2    , -py/2        ,padx,pady),padx,pady);
 576   Info(" 48-  1","sec=%i padx=%3i pady=%3i",p->Loc2Pad(-sx/2-dz , -py/2        ,padx,pady),padx,pady);
 577   Info(" 49-  1","sec=%i padx=%3i pady=%3i",p->Loc2Pad(-sx/2    , -py/2        ,padx,pady),padx,pady);
 578   Info(" 96-  1","sec=%i padx=%3i pady=%3i",p->Loc2Pad( sx/2    , -py/2        ,padx,pady),padx,pady);
 579   Info(" 97-  1","sec=%i padx=%3i pady=%3i",p->Loc2Pad( sx/2+dz , -py/2        ,padx,pady),padx,pady);
 580   Info("144-  1","sec=%i padx=%3i pady=%3i",p->Loc2Pad( px/2    , -py/2        ,padx,pady),padx,pady);
 581   cout<<endl;
 582   Info("  1- 80","sec=%i padx=%3i pady=%3i",p->Loc2Pad(-px/2    , -dz/2        ,padx,pady),padx,pady);
 583   Info(" 48- 80","sec=%i padx=%3i pady=%3i",p->Loc2Pad(-sx/2-dz , -dz/2        ,padx,pady),padx,pady);
 584   Info(" 49- 80","sec=%i padx=%3i pady=%3i",p->Loc2Pad(-sx/2    , -dz/2        ,padx,pady),padx,pady);
 585   Info(" 96- 80","sec=%i padx=%3i pady=%3i",p->Loc2Pad( sx/2    , -dz/2        ,padx,pady),padx,pady);
 586   Info(" 97- 80","sec=%i padx=%3i pady=%3i",p->Loc2Pad( sx/2+dz , -dz/2        ,padx,pady),padx,pady);
 587   Info("144- 80","sec=%i padx=%3i pady=%3i",p->Loc2Pad( px/2    , -dz/2        ,padx,pady),padx,pady);
 588   cout<<endl;
 589   Info("  1- 81","sec=%i padx=%3i pady=%3i",p->Loc2Pad(-px/2    ,  dz/2        ,padx,pady),padx,pady);
 590   Info(" 48- 81","sec=%i padx=%3i pady=%3i",p->Loc2Pad(-sx/2-dz ,  dz/2        ,padx,pady),padx,pady);
 591   Info(" 49- 81","sec=%i padx=%3i pady=%3i",p->Loc2Pad(-sx/2    ,  dz/2        ,padx,pady),padx,pady);
 592   Info(" 96- 81","sec=%i padx=%3i pady=%3i",p->Loc2Pad( sx/2    ,  dz/2        ,padx,pady),padx,pady);
 593   Info(" 97- 81","sec=%i padx=%3i pady=%3i",p->Loc2Pad( sx/2+dz ,  dz/2        ,padx,pady),padx,pady);
 594   Info("144- 81","sec=%i padx=%3i pady=%3i",p->Loc2Pad( px/2    ,  dz/2        ,padx,pady),padx,pady);
 595   cout<<endl;
 596   Info("  1-160","sec=%i padx=%3i pady=%3i",p->Loc2Pad(-px/2    ,  py/2        ,padx,pady),padx,pady);
 597   Info(" 48-160","sec=%i padx=%3i pady=%3i",p->Loc2Pad(-sx/2-dz ,  py/2        ,padx,pady),padx,pady);
 598   Info(" 49-160","sec=%i padx=%3i pady=%3i",p->Loc2Pad(-sx/2    ,  py/2        ,padx,pady),padx,pady);
 599   Info(" 96-160","sec=%i padx=%3i pady=%3i",p->Loc2Pad( sx/2    ,  py/2        ,padx,pady),padx,pady);
 600   Info(" 97-160","sec=%i padx=%3i pady=%3i",p->Loc2Pad( sx/2+dz ,  py/2        ,padx,pady),padx,pady);
 601   Info("144-160","sec=%i padx=%3i pady=%3i",p->Loc2Pad( px/2    ,  py/2        ,padx,pady),padx,pady);
 602   cout<<endl;
 603   Info(" 73-160","sec=%i padx=%3i pady=%3i",p->Loc2Pad(    0    ,  py/2      ,padx,pady),padx,pady);
 604   Info(" 73- 81","sec=%i padx=%3i pady=%3i",p->Loc2Pad(    0    ,  dz/2      ,padx,pady),padx,pady);
 605   Info("0-0dead","sec=%i padx=%3i pady=%3i",p->Loc2Pad(    0    ,   0        ,padx,pady),padx,pady);
 606   Info(" 73- 80","sec=%i padx=%3i pady=%3i",p->Loc2Pad(    0    , -dz/2      ,padx,pady),padx,pady);
 607   Info(" 73-  1","sec=%i padx=%3i pady=%3i",p->Loc2Pad(    0    , -py/2      ,padx,pady),padx,pady);
 608   cout<<endl;
 609   p->Pad2Loc(padx=  1,pady=1,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 610   p->Pad2Loc(padx= 48,pady=1,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 611   p->Pad2Loc(padx= 49,pady=1,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 612   p->Pad2Loc(padx= 96,pady=1,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 613   p->Pad2Loc(padx= 97,pady=1,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 614   p->Pad2Loc(padx=144,pady=1,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 615   cout<<endl;
 616   p->Pad2Loc(padx=  1,pady=80,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 617   p->Pad2Loc(padx= 48,pady=80,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 618   p->Pad2Loc(padx= 49,pady=80,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 619   p->Pad2Loc(padx= 96,pady=80,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 620   p->Pad2Loc(padx= 97,pady=80,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 621   p->Pad2Loc(padx=144,pady=80,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 622   cout<<endl;
 623   p->Pad2Loc(padx=  1,pady=81,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 624   p->Pad2Loc(padx= 48,pady=81,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 625   p->Pad2Loc(padx= 49,pady=81,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 626   p->Pad2Loc(padx= 96,pady=81,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 627   p->Pad2Loc(padx= 97,pady=81,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 628   p->Pad2Loc(padx=144,pady=81,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 629   cout<<endl;
 630   p->Pad2Loc(padx=  1,pady=160,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 631   p->Pad2Loc(padx= 48,pady=160,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 632   p->Pad2Loc(padx= 49,pady=160,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 633   p->Pad2Loc(padx= 96,pady=160,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 634   p->Pad2Loc(padx= 97,pady=160,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 635   p->Pad2Loc(padx=144,pady=160,x,y);  cout<<"padx="<<padx<<" pady="<<pady<<" x="<<x<<" y="<<y<<endl;
 636 }//void TestSeg()
 637 //__________________________________________________________________________________________________
 638 void TestMenu()
 639 {
 640   TControlBar *pMenu = new TControlBar("vertical","RICH test");
 641   pMenu->AddButton("Test segmentation",  "TestSeg()",         "Test AliRICHParam segmentation methods");
 642   pMenu->AddButton("Test transform",     "TestTransform()",   "Test ALiRICHChamber methods");
 643   pMenu->AddButton("Test response",      "TestResponse()",    "Test AliRICHParam response methods");
 644   pMenu->AddButton("Test sdigits",       "TestSD()",          "Create test set of sdigits");
 645   pMenu->AddButton("Test digits OLD",    "TestDigitsOLD()",   "Create test set of OLD digits");
 646   pMenu->AddButton("Test clusters",      "TestC()",           "Create test set of clusters");
 647   pMenu->Show();
 648 }//TestMenu()
 649 //__________________________________________________________________________________________________
 650 void GeoMenu()
 651 {
 652   TControlBar *pMenu = new TControlBar("vertical","RICH draw");
 653   pMenu->AddButton("RICH Isometry", "gMC->Gdraw(\"ALIC\", 60,40,0, 10,10, 0.008,0.008)","Draws ALIC volume in isometry");
 654   pMenu->AddButton("RICH Front XY", "gMC->Gdraw(\"ALIC\", 0,0,0, 10,10, 0.01,0.01)","Draws ALIC volume in XY view");
 655   pMenu->AddButton("RICH Side YZ",  "gMC->Gdraw(\"ALIC\",90,180, 0, 10,10, 0.01,0.01)","Draws ALIC volume in YZ view");
 656   pMenu->AddButton("RICH Top XZ",   "gMC->Gdraw(\"ALIC\",90, 90, 0, 10,10, 0.01,0.01)","Draws ALIC volume in XZ view");
 657   pMenu->AddButton("Module Isometry","gMC->Gdraw(\"SRIC\", 30,60,0, 10,10, 0.1,0.1)","Draws SRIC volume in isometry");
 658   pMenu->AddButton("Module Front XY","gMC->Gdraw(\"SRIC\", 0,0,0, 10,10, 0.1,0.1)","Draws SRIC volume in XY view");
 659   pMenu->AddButton("Module Top XZ", "gMC->Gdraw(\"SRIC\",90, 90, 0, 10,10, 0.1,0.1)","Draws SRIC volume in XZ view");
 660   pMenu->AddButton("ALICE Tree", "((TGeant3*)gMC)->Gdtree(\"ALIC\")","Draws ALICE tree");
 661   pMenu->AddButton("RICH Tree",  "((TGeant3*)gMC)->Gdtree(\"RICH\")","Draws RICH tree");
 662   pMenu->AddButton("Geo test",  "GeoTest()",   "Draw RICH geo as a macro");
 663   pMenu->AddButton("Print ref", "PrintGeo()",  "Print RICH chambers default position");
 664   pMenu->AddButton("AliRICH::Print", "r->Print();", "Print RICH chambers default position");
 665   pMenu->AddButton("Test transform","TestTransform()","Test L2G and G2L methods");
 666   pMenu->Show();
 667 }//GeoMenu()
 668 //__________________________________________________________________________________________________
 669 void menu()
 670 {
 671   TControlBar *pMenu = new TControlBar("vertical","RICH main");
 672
 673   if(ReadAlice()){//it's from file, reconstruct
 674     pMenu->AddButton("hits->sdigits->digits->clusters","MainTrank()","Convert");
 675
 676     pMenu->AddButton("hits->sdigits",    "H_SD()",       "Perform first phase conversion");
 677     pMenu->AddButton("sdigits->digits",  "SD_D()",       "Perform first phase conversion");
 678     pMenu->AddButton("digits->clusters", "D_C()",        "Perform first phase conversion");
 679     pMenu->AddButton("clusters->recos", "C_R()",         "Perform first phase conversion");
 680
 681     pMenu->AddButton("Show","Show()","Shows the structure of events in files");
 682     pMenu->AddButton("OLD specials->sdigits",          "OLD_S_SD()",       "Perform first phase converstion");
 683
 684   }else{//it's aliroot, simulate
 685     pMenu->AddButton("Run",         "a->Run(1)",       "Process!");
 686   }
 687   pMenu->AddButton("Geo submenu",     "GeoMenu()",            "Shows geomentry submenu");
 688   pMenu->AddButton("Geo GUI", "new G3GeometryGUI;","Create instance of G4GeometryGUI");
 689   pMenu->AddButton("Test submenu",    "TestMenu()",            "Shows test submenu");
 690   pMenu->AddButton("Browser",         "new TBrowser;",         "Start ROOT TBrowser");
 691   pMenu->AddButton("Display Fast",    "DisplFast()",           "Display Fast");
 692   pMenu->AddButton("Control Plots",   "ControlPlots()",        "Display some control histograms");
 693   pMenu->AddButton("Debug ON",     "DebugON();",   "Switch debug on-off");
 694   pMenu->AddButton("Debug OFF",    "DebugOFF();",   "Switch debug on-off");
 695   pMenu->AddButton("Quit",            ".q",                    "Close session");
 696   pMenu->Show();
 697   a=gAlice;//for manual manipulation convinience
 698 }//menu()
 699 //__________________________________________________________________________________________________
 700 void DebugOFF(){  Info("DebugOFF","");  a->SetDebug(0);  r->SetDebug(0);  AliLoader::SetDebug(0);}
 701 void DebugON() {  Info("DebugON","");   a->SetDebug(1);  r->SetDebug(1);  AliLoader::SetDebug(1);}