RICH/menu.C

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