RICH/RichMenu.C

   1 AliRun *a;    AliStack *s;  AliRunLoader *al;   TGeoManager *g; //globals for easy manual manipulations
   2 AliRICH   *r; AliLoader    *rl; AliRICHParam *rp;
   3 Bool_t isGeomType=kFALSE;
   4
   5 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
   6 void GetParam()
   7 {
   8   isGeomType=!isGeomType;
   9   if(g) delete g;  if(rp) delete rp; //delete current TGeoManager and AliRICHParam
  10   if(isGeomType) g=TGeoManager::Import("geometry.root");
  11   else           g=TGeoManager::Import("misaligned_geometry.root");
  12   rp=AliRICHParam::Instance();
  13 }
  14 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  15 void RichMenu()
  16 {
  17   TString status="Status: ";
  18   if(gSystem->IsFileInIncludePath("galice.root")){
  19     status+="galice.root found";
  20     al=AliRunLoader::Open();                                                //try to open galice.root from current dir
  21     if(gAlice) delete gAlice;                                               //in case we execute this in aliroot delete default AliRun object
  22     al->LoadgAlice(); a=al->GetAliRun();                                    //take new AliRun object from galice.root
  23     rl=al->GetDetectorLoader("RICH");  r=(AliRICH*)a->GetDetector("RICH");  //get RICH object from galice.root
  24
  25     status+=Form(" with %i event(s)",al->GetNumberOfEvents()); status+=(r)? " with RICH": " without RICH";
  26   }else
  27     status+="No galice.root";
  28
  29   GetParam();
  30
  31   TControlBar *pMenu = new TControlBar("horizontal",status.Data(),0,0);
  32     pMenu->AddButton("                     ","","");
  33     pMenu->AddButton(" General  "           ,"General()"  ,"general items which do not depend on any files");
  34     pMenu->AddButton("                     ",""           ,"");
  35     pMenu->AddButton(" Sim data "           ,"SimData()"  ,"items which expect to have simulated files"    );
  36     pMenu->AddButton("                     ",""           ,"");
  37     pMenu->AddButton(" Raw data "           ,"RawData()"  ,"items which expect to have raw files"          );
  38     pMenu->AddButton("                     ","print()"    ,"");
  39     pMenu->AddButton("Test"                 ,"Test()"     ,"all test utilities");
  40     pMenu->AddButton("                     ","GetParam()" ,"");
  41     pMenu->AddButton("Quit"                 ,".q"         ,"close session"                                 );
  42   pMenu->Show();
  43 }
  44 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  45 void General()
  46 {
  47   TControlBar *pMenu = new TControlBar("vertical","Sim data",60,50);
  48     pMenu->AddButton("Debug ON"           ,"don();"                    ,"Switch debug on-off"                        );
  49     pMenu->AddButton("Debug OFF"          ,"doff();"                   ,"Switch debug on-off"                        );
  50     pMenu->AddButton("Geo GUI"            ,"geo();"                    ,"Shows geometry"                             );
  51     pMenu->AddButton("Browser"            ,"new TBrowser;"             ,"Start ROOT TBrowser"                        );
  52   pMenu->Show();
  53 }//menu()
  54 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  55 void SimData()
  56 {
  57   TControlBar *pMenu = new TControlBar("vertical","Sim",190,50);
  58     pMenu->AddButton("Display ALL chambers"            ,"ed();"     , "Display Fast");
  59     pMenu->AddButton("HITS QA"                         ,"hqa()"     ,"QA plots for hits: hqa()");
  60
  61     pMenu->AddButton("Print hits"                      ,"hp();"      ,"To print hits:     hp(evt)");
  62     pMenu->AddButton("Print sdigits"                   ,"sp();"      ,"To print sdigits:  sp(evt)");
  63     pMenu->AddButton("Print digits"                    ,"dp();"      ,"To print digits:   dp(evt)");
  64     pMenu->AddButton("Print clusters"                  ,"cp();"      ,"To print clusters: cp(evt)");
  65
  66   pMenu->Show();
  67 }
  68 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  69 void RawData()
  70 {
  71   TControlBar *pMenu = new TControlBar("vertical","Raw",350,50);
  72     pMenu->AddButton("ESD print"                       ,"ep();"                  ,"To print ESD info: ep()"         );
  73     pMenu->AddButton("ESD QA"                          ,"eq();"                  ,"To draw ESD hists: eq()"         );
  74     pMenu->AddButton("Clusters print"                  ,"cp();"                  ,"To print clusters: cp()"         );
  75     pMenu->AddButton("Clusters QA"                     ,"cq();"                  ,"To draw clusters hists: cq()"    );
  76     pMenu->AddButton("Print Matrix"                    ,"mp();"                  ,"To print prob matrix: mp()"      );
  77     pMenu->AddButton("Print occupancy"                 ,"r->OccupancyPrint(-1);" ,"To print occupancy"              );
  78     pMenu->AddButton("Print event summary  "           ,"r->SummaryOfEvent();"   ,"To print a summary of the event" );
  79   pMenu->Show();
  80 }//RawData
  81 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  82 void Test()
  83 {
  84   TControlBar *pMenu = new TControlBar("vertical","Test",400,50);
  85     pMenu->AddButton("Hits->Digits"       ,"thd();"                    ,"test hits->sdigits->digits"                 );
  86     pMenu->AddButton("Segmentation"       ,"ts()"                      ,"test segmentation methods"                  );
  87     pMenu->AddButton("Test response"      ,"AliRICHParam::TestResp();" ,"Test AliRICHParam response methods"         );
  88     pMenu->AddButton("Test transformation","AliRICHParam::TestTrans();","Test AliRICHParam transformation methods"   );
  89     pMenu->AddButton("Test Recon"         ,"rec();"                    ,"Test AliRICHRecon"                          );
  90   pMenu->Show();
  91 }//Test()
  92 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
  93
  94
  95 void doff(){  Printf("DebugOFF");  AliLog::SetGlobalDebugLevel(0);}
  96 void don() {  Printf("DebugON");   AliLog::SetGlobalDebugLevel(AliLog::kDebug);}
  97
  98 void geo (                       ) {gGeoManager->GetTopVolume()->Draw("ogl");}
  99
 100 void du  (                       ) {r->Dump         (   );}                //utility display
 101
 102 void hp  (Int_t evt=0            ) {r->HitPrint  (evt);}   //print hits for requested event
 103 void hq  (                       ) {r->HitQA     (   );}   //hits QA plots for all events
 104 void sp  (Int_t evt=0            ) {r->SdiPrint  (evt);}   //print sdigits for requested event
 105 void sq  (Int_t evt=0            ) {r->SdiPrint  (evt);}   //print sdigits for requested event
 106 void dp  (Int_t evt=0            ) {r->DigPrint  (evt);}   //print digits for requested event
 107 void dq  (                       ) {AliRICHReconstructor::DigQA     (al );}   //digits QA plots for all events
 108 void cp  (Int_t evt=0            ) {r->CluPrint  (evt);                   }   //print clusters for requested event
 109 void cq  (                       ) {AliRICHReconstructor::CluQA     (al );}   //clusters QA plots for all events
 110
 111 void ep  (                       ) {AliRICHTracker::EsdQA(1);            }
 112 void eq  (                       ) {AliRICHTracker::EsdQA();             }
 113 void mp  (Double_t probCut=0.7   ) {AliRICHTracker::MatrixPrint(probCut);}
 114
 115
 116 void t   (Int_t evt=0          )   {AliRICHParam::Stack(evt);}
 117 void tid (Int_t tid,Int_t evt=0)   {AliRICHParam::Stack(evt,tid);}
 118
 119
 120 Int_t nem (Int_t evt=0) {AliRICHParam::StackCount(kElectron  ,evt);} //utility number of electrons
 121 Int_t nep (Int_t evt=0) {AliRICHParam::StackCount(kPositron  ,evt);} //utility number of positrons
 122 Int_t nmup(Int_t evt=0) {AliRICHParam::StackCount(kMuonPlus  ,evt);} //utility number of positive muons
 123 Int_t nmum(Int_t evt=0) {AliRICHParam::StackCount(kMuonMinus ,evt);} //utility number of negative muons
 124 Int_t npi0(Int_t evt=0) {AliRICHParam::StackCount(kPi0       ,evt);} //utility number of neutral pions
 125 Int_t npip(Int_t evt=0) {AliRICHParam::StackCount(kPiPlus    ,evt);} //utility number of positive pions
 126 Int_t npim(Int_t evt=0) {AliRICHParam::StackCount(kPiMinus   ,evt);} //utility number of negative pions
 127 Int_t nk0 (Int_t evt=0) {AliRICHParam::StackCount(kK0        ,evt);} //utility number of neutral kaons
 128 Int_t nkp (Int_t evt=0) {AliRICHParam::StackCount(kKPlus     ,evt);} //utility number of positive kaons
 129 Int_t nkm (Int_t evt=0) {AliRICHParam::StackCount(kKMinus    ,evt);} //utility number of negative kaons
 130 Int_t npp (Int_t evt=0) {AliRICHParam::StackCount(kProton    ,evt);} //utility number of protons
 131 Int_t npm (Int_t evt=0) {AliRICHParam::StackCount(kProtonBar ,evt);} //utility number of antiprotons
 132 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 133
 134 void tst()
 135 {
 136   Int_t ch=3;
 137   TClonesArray hitLst("AliRICHHit");
 138   TClonesArray sdiLst("AliRICHDigit");
 139   TObjArray    digLst;                  for(Int_t i=0;i<7;i++) digLst.AddAt(new TClonesArray("AliRICHDigit"),i);
 140   TClonesArray cluLst("AliRICHCluster");
 141
 142   Int_t iHitCnt=0;  new(hitLst[iHitCnt++]) AliRICHHit(ch,200e-9,kProton,1, 8.40 , 60.14); //c,e,pid,tid,xl,yl
 143 //                    new(hitLst[iHitCnt++]) AliRICHHit(ch,e*1e-9,kProton,2,x,y); //c,e,pid,tid,xl,yl,x,y,z
 144
 145                                                                      Printf("HIT------HIT---------HIT--------HIT------HIT------HIT");hitLst.Print();Printf("");
 146   AliRICHv1::Hit2Sdi(&hitLst,&sdiLst);                               Printf("SDI------DIG---------SDI--------SDI------SDI------SDI");sdiLst.Print();Printf("");
 147   AliRICHDigitizer::Sdi2Dig(&sdiLst,&digLst);                        Printf("DIG------DIG---------DIG--------DIG------DIG------DIG");digLst.Print();Printf("");
 148   AliRICHReconstructor::Dig2Clu((TClonesArray*)digLst[ch],&cluLst);  Printf("CLU------CLU---------CLU--------CLU------CLU------CLU");cluLst.Print();Printf("");
 149
 150
 151 }
 152
 153
 154 void print()
 155 {
 156
 157   Double_t r2d=TMath::RadToDeg();
 158
 159   Double_t x=AliRICHDigit::SizeAllX(),y=AliRICHDigit::SizeAllY();
 160
 161   TVector3 c6lt=rp->Lors2Mars(6,0,y);                                              TVector3 c6rt=rp->Lors2Mars(6,x,y);
 162                                        TVector3 c6ce=rp->Lors2Mars(6,x/2,y/2);
 163   TVector3 c6lb=rp->Lors2Mars(6,0,0);                                              TVector3 c6rb=rp->Lors2Mars(6,x,0);
 164 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 165   TVector3 c5lt=rp->Lors2Mars(5,0,y);                                              TVector3 c5rt=rp->Lors2Mars(5,x,y);
 166                                        TVector3 c5ce=rp->Lors2Mars(5,x/2,y/2);
 167   TVector3 c5lb=rp->Lors2Mars(5,0,0);                                              TVector3 c5rb=rp->Lors2Mars(5,x,0);
 168 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 169   TVector3 c4lt=rp->Lors2Mars(4,0,y);                                              TVector3 c4rt=rp->Lors2Mars(4,x,y);
 170                                        TVector3 c4ce=rp->Lors2Mars(4,x/2,y/2);
 171   TVector3 c4lb=rp->Lors2Mars(4,0,0);                                              TVector3 c4rb=rp->Lors2Mars(4,x,0);
 172 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 173   TVector3 c3lt=rp->Lors2Mars(3,0,y);                                              TVector3 c3rt=rp->Lors2Mars(3,x,y);
 174                                        TVector3 c3ce=rp->Lors2Mars(3,x/2,y/2);
 175   TVector3 c3lb=rp->Lors2Mars(3,0,0);                                              TVector3 c3rb=rp->Lors2Mars(3,x,0);
 176 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 177   TVector3 c2lt=rp->Lors2Mars(2,0,y);                                              TVector3 c2rt=rp->Lors2Mars(2,x,y);
 178                                        TVector3 c2ce=rp->Lors2Mars(2,x/2,y/2);
 179   TVector3 c2lb=rp->Lors2Mars(2,0,0);                                              TVector3 c2rb=rp->Lors2Mars(2,x,0);
 180 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 181   TVector3 c1lt=rp->Lors2Mars(1,0,y);                                              TVector3 c1rt=rp->Lors2Mars(1,x,y);
 182                                        TVector3 c1ce=rp->Lors2Mars(1,x/2,y/2);
 183   TVector3 c1lb=rp->Lors2Mars(1,0,0);                                              TVector3 c1rb=rp->Lors2Mars(1,x,0);
 184 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 185   TVector3 c0lt=rp->Lors2Mars(0,0,y);                                              TVector3 c0rt=rp->Lors2Mars(0,x,y);
 186                                        TVector3 c0ce=rp->Lors2Mars(0,x/2,y/2);
 187   TVector3 c0lb=rp->Lors2Mars(0,0,0);                                              TVector3 c0rb=rp->Lors2Mars(0,x,0);
 188
 189
 190   Printf("\n\n\n");
 191
 192   Printf("_______________________________   _______________________________");
 193   Printf("|%7.2f               %7.2f|   |%7.2f               %7.2f|",c6lt.Mag()          ,c6rt.Mag()        ,c5lt.Mag()        ,c5rt.Mag()                 );
 194   Printf("|%7.2f               %7.2f|   |%7.2f               %7.2f|",c6lt.Theta()*r2d    ,c6rt.Theta()*r2d  ,c5lt.Theta()*r2d  ,c5rt.Theta()*r2d           );
 195   Printf("|%7.2f               %7.2f|   |%7.2f               %7.2f|",c6lt.Phi()*r2d      ,c6rt.Phi()*r2d    ,c5lt.Phi()*r2d    ,c5rt.Phi()*r2d             );
 196   Printf("|                             |   |                             |");
 197   Printf("|         %7.2f             |   |         %7.2f             | Sensitive area  (130.60,126.16)"  ,c6ce.Mag()           ,c5ce.Mag()                );
 198   Printf("|         %7.2f             |   |         %7.2f             | Lors Center     ( 65.30, 63.08)"  ,c6ce.Theta()*r2d     ,c5ce.Theta()*r2d          );
 199   Printf("|         %7.2f             |   |         %7.2f             |"                                  ,c6ce.Phi()*r2d       ,c5ce.Phi()*r2d            );
 200   Printf("|                             |   |                             |");
 201   Printf("|%7.2f               %7.2f|   |%7.2f               %7.2f|",c6lb.Mag()          ,c6rb.Mag()        ,c5lb.Mag()        ,c5rb.Mag()                 );
 202   Printf("|%7.2f               %7.2f|   |%7.2f               %7.2f|",c6lb.Theta()*r2d    ,c6rb.Theta()*r2d  ,c5lb.Theta()*r2d  ,c5rb.Theta()*r2d           );
 203   Printf("|%7.2f               %7.2f|   |%7.2f               %7.2f|",c6lb.Phi()*r2d      ,c6rb.Phi()*r2d    ,c5lb.Phi()*r2d    ,c5rb.Phi()*r2d             );
 204   Printf("-------------------------------   -------------------------------");
 205   Printf("");
 206   Printf("_______________________________   _______________________________   _______________________________");
 207   Printf("|%7.2f               %7.2f|   |%7.2f               %7.2f|   |%7.2f               %7.2f|",c4lt.Mag()      ,c4rt.Mag()      ,c3lt.Mag()      ,c3rt.Mag()      ,c2lt.Mag()       ,c2rt.Mag()      );
 208   Printf("|%7.2f               %7.2f|   |%7.2f               %7.2f|   |%7.2f               %7.2f|",c4lt.Theta()*r2d,c4rt.Theta()*r2d,c3lt.Theta()*r2d,c3rt.Theta()*r2d,c2lt.Theta()*r2d ,c2rt.Theta()*r2d);
 209   Printf("|%7.2f               %7.2f|   |%7.2f               %7.2f|   |%7.2f               %7.2f|",c4lt.Phi()*r2d  ,c4rt.Phi()*r2d  ,c3lt.Phi()*r2d  ,c3rt.Phi()*r2d  ,c2lt.Phi()*r2d   ,c2rt.Phi()*r2d  );
 210   Printf("|                             |   |                             |   |                             |");
 211   Printf("|         %7.2f             |   |         %7.2f             |   |         %7.2f             |"     ,c4ce.Mag()      ,c3ce.Mag()      ,c2ce.Mag());
 212   Printf("|         %7.2f             |   |         %7.2f             |   |         %7.2f             |"     ,c4ce.Theta()*r2d,c3ce.Theta()*r2d,c2ce.Theta()*r2d);
 213   Printf("|         %7.2f             |   |         %7.2f             |   |         %7.2f             |"     ,c4ce.Phi()*r2d  ,c3ce.Phi()*r2d  ,c3ce.Phi()*r2d);
 214   Printf("|                             |   |                             |   |                             |");
 215   Printf("|%7.2f               %7.2f|   |%7.2f               %7.2f|   |%7.2f               %7.2f|",c4lb.Mag()      ,c4rb.Mag()      ,c3lb.Mag()      ,c3rb.Mag()      ,c2lt.Mag()       ,c2rt.Mag()      );
 216   Printf("|%7.2f               %7.2f|   |%7.2f               %7.2f|   |%7.2f               %7.2f|",c4lb.Theta()*r2d,c4rb.Theta()*r2d,c3lb.Theta()*r2d,c3rb.Theta()*r2d,c2lt.Theta()*r2d ,c2rt.Theta()*r2d);
 217   Printf("|%7.2f               %7.2f|   |%7.2f               %7.2f|   |%7.2f               %7.2f|",c4lb.Phi()*r2d  ,c4rb.Phi()*r2d  ,c3lb.Phi()*r2d  ,c3rb.Phi()*r2d  ,c2lt.Phi()*r2d   ,c2rt.Phi()*r2d  );
 218   Printf("-------------------------------   -------------------------------   -------------------------------");
 219   Printf("");
 220   Printf("                                  _______________________________   _______________________________");
 221   Printf("                                  |%7.2f               %7.2f|   |%7.2f               %7.2f|",c1lt.Mag()      ,c1rt.Mag()      ,c0lt.Mag()      ,c0rt.Mag()      );
 222   Printf("                                  |%7.2f               %7.2f|   |%7.2f               %7.2f|",c1lt.Theta()*r2d,c1rt.Theta()*r2d,c0lt.Theta()*r2d,c0rt.Theta()*r2d);
 223   Printf("                                  |%7.2f               %7.2f|   |%7.2f               %7.2f|",c1lt.Phi()*r2d  ,c1rt.Phi()*r2d  ,c0lt.Phi()*r2d  ,c0rt.Phi()*r2d  );
 224   Printf("                                  |                             |   |                             |");
 225   Printf("                                  |         %7.2f             |   |         %7.2f             |"     ,c1ce.Mag()      ,c0ce.Mag()      );
 226   Printf("                                  |         %7.2f             |   |         %7.2f             |"     ,c1ce.Theta()*r2d,c0ce.Theta()*r2d);
 227   Printf("                                  |         %7.2f             |   |         %7.2f             |"     ,c1ce.Phi()*r2d  ,c0ce.Phi()*r2d  );
 228   Printf("                                  |                             |   |                             |");
 229   Printf("                                  |%7.2f               %7.2f|   |%7.2f               %7.2f|",c1lb.Mag()      ,c1rb.Mag()      ,c0lb.Mag()      ,c0rb.Mag()       );
 230   Printf("                                  |%7.2f               %7.2f|   |%7.2f               %7.2f|",c1lb.Theta()*r2d,c1rb.Theta()*r2d,c0lb.Theta()*r2d,c0rb.Theta()*r2d );
 231   Printf("                                  |%7.2f               %7.2f|   |%7.2f               %7.2f|",c1lb.Phi()*r2d  ,c1rb.Phi()*r2d  ,c0lb.Phi()*r2d  ,c0rb.Phi()*r2d   );
 232   Printf("                                  -------------------------------   -------------------------------");
 233
 234
 235   Double_t m[3];
 236   for(int i=0;i<1000;i++){
 237     Float_t xout=0,xin=gRandom->Rndm()*130.60;
 238     Float_t yout=0,yin=gRandom->Rndm()*126.16;
 239     Int_t   c=gRandom->Rndm()*6;
 240     rp->Lors2Mars(c,xin,yin,m);
 241     rp->Mars2Lors(c,m,xout,yout);
 242     if( (xin-xout) != 0) Printf("Problem in X");
 243     if( (yin-yout) != 0) Printf("Problem in Y");
 244   }
 245 }//print()
 246
 247
 248 void rec()
 249 {
 250   Double_t th=8*TMath::DegToRad();//gRandom->Rndm()*TMath::PiOver4();
 251   Double_t ph=gRandom->Rndm()*TMath::TwoPi(); 172.51*TMath::DegToRad();
 252   Double_t  x=gRandom->Rndm()*AliRICHDigit::SizeAllX(); //101.59;
 253   Double_t  y=gRandom->Rndm()*2*AliRICHDigit::SizePcY();//38.06;
 254
 255
 256   Double_t ckovMax=TMath::ACos(1./AliRICHRecon::fkRadIdx),ckovSim;
 257   Int_t nSim=0;
 258   while(nSim<3){
 259     ckovSim=gRandom->Rndm()*ckovMax;//0.6468;
 260     nSim=20*TMath::Power(TMath::Sin(ckovSim)/TMath::Sin(ckovMax),2); //scale number of photons
 261   }
 262
 263
 264   TClonesArray *pCluLst=new TClonesArray("AliRICHCluster");
 265   TPolyMarker  *pMipMap=new TPolyMarker();   pMipMap->SetMarkerStyle(8);  pMipMap->SetMarkerColor(kRed);
 266   TPolyMarker  *pPhoMap=new TPolyMarker();   pPhoMap->SetMarkerStyle(4);  pPhoMap->SetMarkerColor(kRed);
 267   TPolyMarker  *pBkgMap=new TPolyMarker();   pBkgMap->SetMarkerStyle(25); pBkgMap->SetMarkerColor(kRed);
 268   TPolyLine    *pRing  =new TPolyLine;                                    pRing->SetLineColor(kGreen);
 269   TPolyLine    *pOld   =new TPolyLine;                                    pOld->SetLineColor(kBlue);
 270
 271   Int_t iCluCnt=0; pMipMap->SetPoint(iCluCnt,x,y); new((*pCluLst)[iCluCnt++])   AliRICHCluster(1,x,y,200); //add mip cluster
 272
 273 //  for(int j=0;j<30;j++){                                                                                   //add bkg photons
 274 //    Float_t bkgX=gRandom->Rndm()*AliRICHDigit::SizeAllX();
 275 //    Float_t bkgY=gRandom->Rndm()*AliRICHDigit::SizeAllY();
 276 //    pBkgMap->SetPoint(iCluCnt,bkgX,bkgY); new((*pCluLst)[iCluCnt++]) AliRICHCluster(1,bkgX,bkgY,35);
 277 //  }
 278
 279
 280
 281
 282   AliRICHRecon    rec; rec.SetTrack(th,ph,x,y);
 283   AliRICHReconOld old; old.SetTrack(th,ph,x,y);
 284
 285   TVector2 pos;
 286   for(int i=0;i<nSim;i++){
 287     rec.TracePhoton(ckovSim,gRandom->Rndm()*2*TMath::Pi(),pos);                                   //add photons
 288     if(AliRICHDigit::IsInDead(pos.X(),pos.Y())) continue;
 289     pPhoMap->SetPoint(iCluCnt,pos.X(),pos.Y());    new((*pCluLst)[iCluCnt++]) AliRICHCluster(1,pos.X(),pos.Y(),35);
 290   }
 291
 292   for(int i=0;i<pCluLst->GetEntries();i++){
 293     AliRICHCluster *pClu=(AliRICHCluster*)pCluLst->At(i);
 294     Printf("diff phi %f ckov new  %f ckov old %f",rec.FindPhotPhi(pClu->X(),pClu->Y())-old.FindPhotPhi(pClu->X(),pClu->Y()),rec.FindPhotCkov(pClu->X(),pClu->Y()),old.FindPhotCkov(pClu->X(),pClu->Y()));
 295   }
 296
 297   Int_t nRec=0,nOld=0;
 298   Double_t ckovRec=rec.CkovAngle(pCluLst,nRec); Double_t err=TMath::Sqrt(rec.CkovSigma2());
 299   Double_t ckovOld=old.CkovAngle(pCluLst,nOld);
 300
 301   Printf("---------------- Now reconstructed --------------------");
 302
 303
 304   for(int j=0;j<100;j++){rec.TracePhoton(ckovRec,j*0.0628,pos); pRing->SetPoint(j,pos.X(),pos.Y());}
 305   for(int j=0;j<100;j++){rec.TracePhoton(ckovOld,j*0.0628,pos); pOld->SetPoint(j,pos.X(),pos.Y());}
 306
 307   new TCanvas;  AliRICHDigit::DrawPc();  pMipMap->Draw(); pPhoMap->Draw(); pBkgMap->Draw(); pRing->Draw();  pOld->Draw();
 308
 309   TLatex txt; txt.SetTextSize(0.03);
 310   txt.DrawLatex(65,127,Form("#theta=%.4f#pm%.5f with %2i #check{C}"                             ,ckovSim, 0.,nSim             ));
 311   txt.DrawLatex(65,122,Form("#theta=%.4f#pm%.5f with %2i #check{C} Old=%.4f with %i #check{C}"  ,ckovRec,err,nRec,ckovOld,nOld));
 312   txt.DrawLatex(0 ,127,Form("#theta=%.2f#circ   #phi=%.2f#circ @(%.2f,%.2f) ",th*TMath::RadToDeg(),ph*TMath::RadToDeg(),x,y));
 313
 314 //  for(int i=0;i<35;i++){
 315 //    Double_t ckov=0.1+i*0.02;
 316 //    Printf("Ckov=%.2f Old=%.3f New=%.3f",ckov,old.FindRingArea(ckov),rec.FindRingArea(ckov));
 317 //  }
 318
 319
 320 }//rec()
 321
 322 AliRICHRecon rec; AliRICHReconOld old;
 323
 324 void aaa()
 325 {
 326   rec.SetTrack(10*TMath::DegToRad(),1,30,60); old.SetTrack(10*TMath::DegToRad(),1,30,60);
 327   TVector2 pos;
 328   Double_t ckovSim=0.6234;
 329   rec.TracePhoton(ckovSim,1,pos);
 330   Printf("ckovSim %f",ckovSim);
 331   double ckovRec=rec.FindPhotCkov(pos.X(),pos.Y());
 332   double ckovOld=old.FindPhotCkov(pos.X(),pos.Y());
 333   Printf("new %f old %f",ckovRec,ckovOld);
 334 }
 335
 336
 337
 338
 339
 340 void ed(Int_t iEvtFrom=-1,Int_t iEvtTo=-1)
 341 {
 342 // Display digits, reconstructed tracks intersections and RICH rings if available
 343
 344
 345   Int_t iNevt=al->GetNumberOfEvents();
 346   if(iEvtFrom< 0    ){iEvtFrom=0;iEvtTo=iNevt-1;}
 347   if(iEvtTo  >=iNevt){           iEvtTo=iNevt-1;}
 348
 349   rl->LoadDigits();
 350
 351   TFile *pEsdFl=TFile::Open("AliESDs.root");     if(!pEsdFl || !pFile->IsOpen()) return;//open AliESDs.root
 352   TTree *pEsdTr=(TTree*) pFile->Get("esdTree");  if(!pEsdTr)                     return;//get ESD tree
 353
 354   AliESD *pEsd=new AliESD;  pTree->SetBranchAddress("ESD", &pEsd);
 355
 356
 357
 358
 359   TPolyMarker  *pDigMap[7]; //digits map
 360   TPolyMarker  *pTrkMap[7]; //TRKxPC intersection map
 361
 362   for(Int_t i=0;i<7;i++){
 363     pDigMap[i]=new TPolyMarker();       pDigMap[i]->SetMarkerStyle(25);   pDigMap[i]->SetMarkerSize(0.5); pDigMap[i]->SetMarkerColor(kGreen);
 364     pTrkMap[i]=new TPolyMarker();       pTrkMap[i]->SetMarkerStyle(4);    pTrkMap[i]->SetMarkerSize(0.5); pTrkMap[i]->SetMarkerColor(kRed);
 365   }
 366
 367   TLatex t;
 368   TCanvas *pC = new TCanvas("RICHDisplay","RICH Display",0,0,1226,900);  pC->Divide(3,3);
 369
 370   for(Int_t iEvt=iEvtFrom;iEvt<=iEvtTo;iEvt++) {                //events loop
 371     pC->cd(3);  t.DrawText(0.2,0.4,Form("Event %i",iEvt));        //print current event number
 372
 373     pEsdTr->GetEntry(iEvt);                              //get ESD for this event
 374     for(Int_t iTrk=0;iTrk<pEsd->GetNumberOfTracks();iTrk++){//ESD tracks loop
 375       AliESDtrack *pTrk = pEsd->GetTrack(iTrk);             //
 376       Float_t th,ph,x,y;
 377       pTrk->GetRICHtrk(th,ph,x,y);
 378
 379     }//ESD tracks loop
 380
 381     al->GetEvent(iEvt);   rl->TreeD()->GetEntry(0); //get digits list
 382     for(Int_t iCh=0;iCh<7;iCh++) {//chambers loop
 383       for(Int_t iDig=0;iDig < r->DigLst(iCh)->GetEntries();iDig++) {      //digits loop
 384         AliRICHDigit *pDig = (AliRICHDigit*)r->DigLst(iCh)->At(iDig);
 385         pDigMap[iCh]->SetPoint(iDig,pDig->LorsX(),pDig->LorsY());
 386       }                                                             //digits loop
 387
 388
 389       if(iCh==6) pC->cd(1); if(iCh==5) pC->cd(2);
 390       if(iCh==4) pC->cd(4); if(iCh==3) pC->cd(5); if(iCh==2) pC->cd(6);
 391                             if(iCh==1) pC->cd(8); if(iCh==0) pC->cd(9);
 392
 393       AliRICHDigit::DrawPc();  pDigMap[iCh]->Draw();
 394     }//chambers loop
 395     pC->Update();
 396     pC->Modified();
 397
 398     if(iEvt<iEvtTo) {gPad->WaitPrimitive();pC->Clear();}
 399
 400
 401
 402   }//events loop
 403   delete pEsf;  pEsdFl->Close();//close AliESDs.root
 404   rl->UnloadDigits();
 405 }//Display()
 406
 407
 408 void AliRICH::HitQA(Double_t cut,Double_t cutele,Double_t cutR)
 409 {
 410 // Provides a set of control plots intended primarily for charged particle flux analisys
 411 // Arguments: cut (GeV)    - cut on momentum of any charged particles but electrons,
 412 //            cetele (GeV) - the same for electrons-positrons
 413 //            cutR (cm)    - cut on production vertex radius (cylindrical system)
 414   gBenchmark->Start("HitsAna");
 415
 416   Double_t cutPantiproton    =cut;
 417   Double_t cutPkaonminus     =cut;
 418   Double_t cutPpionminus     =cut;
 419   Double_t cutPmuonminus     =cut;
 420   Double_t cutPpositron      =cutele;
 421
 422   Double_t cutPelectron      =cutele;
 423   Double_t cutPmuonplus      =cut;
 424   Double_t cutPpionplus      =cut;
 425   Double_t cutPkaonplus      =cut;
 426   Double_t cutPproton        =cut;
 427
 428
 429   TH2F *pEleHitRZ    =new TH2F("EleHitRZ"    ,Form("e^{+} e^{-} hit %s;z[cm];R[cm]" ,GetName())     , 400,-300,300 ,400,-500,500);   //R-z plot 0cm<R<550cm -300cm<z<300cm
 430   TH2F *pEleHitRP    =new TH2F("EleHitRP"    ,Form("e^{+} e^{-} hit %s;p[GeV];R[cm]",GetName())     ,1000,-1  ,1   ,400,   0,550);   //R-p plot 0cm<R<550cm -1GeV<p<1GeV
 431   TH1F *pEleAllP     =new TH1F("EleAllP"     ,     "e^{+} e^{-} all;p[GeV]"                         ,1000,-1  ,1                );
 432   TH1F *pEleHitP     =new TH1F("EleHitP"     ,Form("e^{+} e^{-} hit %s;p[GeV]"      ,GetName())     ,1000,-1  ,1                );
 433   TH1F *pMuoHitP     =new TH1F("MuoHitP"     ,Form("#mu^{-} #mu^{+} hit %s;p[GeV]"  ,GetName())     ,1000,-4  ,4                );
 434   TH1F *pPioHitP     =new TH1F("PioHitP"     ,Form("#pi^{-} #pi^{+} hit %s;p[GeV]"  ,GetName())     ,1000,-4  ,4                );
 435   TH1F *pKaoHitP     =new TH1F("KaoHitP"     ,Form("K^{-} K^{+} hit %s;p[GeV]"      ,GetName())     ,1000,-4  ,4                );
 436   TH1F *pProHitP     =new TH1F("ProHitP"     ,Form("p^{-} p^{+} hit %s;p[GeV]"      ,GetName())     ,1000,-4  ,4                );
 437   TH2F *pFlux        =new TH2F("flux"        ,Form("%s flux with Rvertex<%.1fcm"    ,GetName(),cutR),10  ,-5  ,5   , 10,0    ,10); //special text hist
 438   TH2F *pVertex      =new TH2F("vertex"      ,Form("%s 2D vertex of RICH hit;x;y"   ,GetName())     ,120 ,0   ,600 ,120,0    ,600); //special text hist
 439   TH1F *pRho         =new TH1F("rho"         ,Form("%s r of RICH hit"               ,GetName())     ,600 ,0   ,600); //special text hist
 440   pFlux->SetStats(0);
 441   pFlux->GetXaxis()->SetBinLabel(1 ,Form("p^{-}>%.3fGeV/c"   ,cutPantiproton));
 442   pFlux->GetXaxis()->SetBinLabel(2 ,Form("K^{-}>%.3fGeV/c"   ,cutPkaonminus ));
 443   pFlux->GetXaxis()->SetBinLabel(3 ,Form("#pi^{-}>%.3fGeV/c" ,cutPpionminus ));
 444   pFlux->GetXaxis()->SetBinLabel(4 ,Form("#mu^{-}>%.3fGeV/c" ,cutPmuonminus ));
 445   pFlux->GetXaxis()->SetBinLabel(5 ,Form("e^{+}>%.3fGeV/c"   ,cutPpositron  ));
 446
 447   pFlux->GetXaxis()->SetBinLabel(6 ,Form("e^{-}>%.3fGeV/c"   ,cutPelectron  ));
 448   pFlux->GetXaxis()->SetBinLabel(7 ,Form("#mu^{+}>%.3fGeV/c" ,cutPmuonplus  ));
 449   pFlux->GetXaxis()->SetBinLabel(8 ,Form("#pi^{+}>%.3fGeV/c" ,cutPpionplus  ));
 450   pFlux->GetXaxis()->SetBinLabel(9 ,Form("K^{+}>%.3fGeV/c"   ,cutPkaonplus  ));
 451   pFlux->GetXaxis()->SetBinLabel(10,Form("p^{+}>%.3fGeV/c"   ,cutPproton    ));
 452
 453   pFlux->GetYaxis()->SetBinLabel(1,"sum");
 454   pFlux->GetYaxis()->SetBinLabel(2,"ch1");
 455   pFlux->GetYaxis()->SetBinLabel(3,"ch2");
 456   pFlux->GetYaxis()->SetBinLabel(4,"ch3");
 457   pFlux->GetYaxis()->SetBinLabel(5,"ch4");
 458   pFlux->GetYaxis()->SetBinLabel(6,"ch5");
 459   pFlux->GetYaxis()->SetBinLabel(7,"ch6");
 460   pFlux->GetYaxis()->SetBinLabel(8,"ch7");
 461   pFlux->GetYaxis()->SetBinLabel(9,"prim");
 462   pFlux->GetYaxis()->SetBinLabel(10,"tot");
 463
 464 //end of hists definition
 465
 466   Int_t iNevents=fLoader->GetRunLoader()->GetAliRun()->GetEventsPerRun(),iCntPrimParts=0,iCntTotParts=0;
 467 //load all needed trees
 468   fLoader->LoadHits();
 469   fLoader->GetRunLoader()->LoadHeader();
 470   fLoader->GetRunLoader()->LoadKinematics();
 471
 472   for(Int_t iEvtN=0;iEvtN < iNevents;iEvtN++){//events loop
 473     fLoader->GetRunLoader()->GetEvent(iEvtN);
 474     AliInfo(Form(" %i event processes",fLoader->GetRunLoader()->GetEventNumber()));
 475     AliStack *pStack= fLoader->GetRunLoader()->Stack();
 476
 477     for(Int_t iParticleN=0;iParticleN<pStack->GetNtrack();iParticleN++){//stack loop
 478       TParticle *pPart=pStack->Particle(iParticleN);
 479
 480       if(iParticleN%10000==0) AliInfo(Form(" %i particles read",iParticleN));
 481
 482       switch(pPart->GetPdgCode()){
 483         case kProtonBar: pFlux->Fill(-4.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-4.5,8); break;
 484         case kKMinus:    pFlux->Fill(-3.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-3.5,8); break;
 485         case kPiMinus:   pFlux->Fill(-2.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-2.5,8); break;
 486         case kMuonMinus: pFlux->Fill(-1.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-1.5,8); break;
 487         case kPositron:  pFlux->Fill(-0.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-0.5,8); pEleAllP->Fill(-pPart->P()); break;
 488
 489         case kElectron:  pFlux->Fill( 0.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 0.5,8); pEleAllP->Fill( pPart->P()); break;
 490         case kMuonPlus:  pFlux->Fill( 1.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 1.5,8); break;
 491         case kPiPlus:    pFlux->Fill( 2.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 2.5,8); break;
 492         case kKPlus:     pFlux->Fill( 3.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 3.5,8); break;
 493         case kProton:    pFlux->Fill( 4.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 4.5,8); break;
 494       }//switch
 495     }//stack loop
 496 //now hits analiser
 497     for(Int_t iEntryN=0;iEntryN < fLoader->TreeH()->GetEntries();iEntryN++){//TreeH loop
 498       fLoader->TreeH()->GetEntry(iEntryN);                                  //get current entry (prim)
 499       for(Int_t iHitN=0;iHitN < Hits()->GetEntries();iHitN++){//hits loop
 500         AliRICHHit *pHit = (AliRICHHit*)Hits()->At(iHitN);            //get current hit
 501         TParticle  *pPart=pStack->Particle(pHit->GetTrack());      //get stack particle which produced the current hit
 502
 503         if(pPart->GetPDG()->Charge()!=0&&pPart->Rho()>0.1) pVertex->Fill(pPart->Vx(),pPart->Vy()); //safe margin for sec.
 504         if(pPart->GetPDG()->Charge()!=0) pRho->Fill(pPart->Rho()); //safe margin for sec.
 505         if(pPart->R()>cutR) continue;                                   //cut on production radius (cylindrical system)
 506
 507         switch(pPart->GetPdgCode()){
 508           case kProtonBar: if(pPart->P()>cutPantiproton) {pProHitP->Fill(-pPart->P()); pFlux->Fill(-4.5,pHit->Ch());}break;
 509           case kKMinus   : if(pPart->P()>cutPkaonminus)  {pKaoHitP->Fill(-pPart->P()); pFlux->Fill(-3.5,pHit->Ch());}break;
 510           case kPiMinus  : if(pPart->P()>cutPpionminus)  {pPioHitP->Fill(-pPart->P()); pFlux->Fill(-2.5,pHit->Ch());}break;
 511           case kMuonMinus: if(pPart->P()>cutPmuonminus)  {pMuoHitP->Fill(-pPart->P()); pFlux->Fill(-1.5,pHit->Ch());}break;
 512           case kPositron : if(pPart->P()>cutPpositron)   {pEleHitP->Fill(-pPart->P()); pFlux->Fill(-0.5,pHit->Ch());
 513                pEleHitRP->Fill(-pPart->P(),pPart->R());  pEleHitRZ->Fill(pPart->Vz(),pPart->R()); }break;
 514
 515           case kElectron : if(pPart->P()>cutPelectron)   {pEleHitP->Fill( pPart->P()); pFlux->Fill( 0.5,pHit->Ch());
 516                pEleHitRP->Fill( pPart->P(),pPart->R());  pEleHitRZ->Fill(pPart->Vz(),pPart->R()); }break;
 517           case kMuonPlus : if(pPart->P()>cutPmuonplus)   {pMuoHitP->Fill( pPart->P()); pFlux->Fill( 1.5,pHit->Ch());}break;
 518           case kPiPlus   : if(pPart->P()>cutPpionplus)   {pPioHitP->Fill( pPart->P()); pFlux->Fill( 2.5,pHit->Ch());}break;
 519           case kKPlus    : if(pPart->P()>cutPkaonplus)   {pKaoHitP->Fill( pPart->P()); pFlux->Fill( 3.5,pHit->Ch());}break;
 520           case kProton   : if(pPart->P()>cutPproton)     {pProHitP->Fill( pPart->P()); pFlux->Fill( 4.5,pHit->Ch());}break;
 521         }
 522       }//hits loop
 523     }//TreeH loop
 524     iCntPrimParts +=pStack->GetNprimary();
 525     iCntTotParts  +=pStack->GetNtrack();
 526   }//events loop
 527 //unload all loaded staff
 528   fLoader->UnloadHits();
 529   fLoader->GetRunLoader()->UnloadHeader();
 530   fLoader->GetRunLoader()->UnloadKinematics();
 531 //Calculater some sums
 532   Stat_t sum=0;
 533 //sum row, sum over rows
 534   for(Int_t i=1;i<=pFlux->GetNbinsX();i++){
 535     sum=0; for(Int_t j=2;j<=8;j++)    sum+=pFlux->GetBinContent(i,j);
 536     pFlux->SetBinContent(i,1,sum);
 537   }
 538
 539 //display everything
 540   new TCanvas("canvas1",Form("Events %i Nprims=%i Nparticles=%i",iNevents,iCntPrimParts,iCntTotParts),1000,900); pFlux->Draw("text");  gPad->SetGrid();
 541 //total prims and particles
 542   TLatex latex; latex.SetTextSize(0.02);
 543   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));
 544   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));
 545   for(Int_t iCh=0;iCh<7;iCh++) {
 546     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));
 547   }
 548   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));
 549
 550   new TCanvas("cEleAllP"   ,"e" ,200,100); pEleAllP->Draw();
 551   new TCanvas("cEleHitRP"  ,"e" ,200,100); pEleHitRP->Draw();
 552   new TCanvas("cEleHitRZ"  ,"e" ,200,100); pEleHitRZ->Draw();
 553   new TCanvas("cEleHitP"   ,"e" ,200,100); pEleHitP->Draw();
 554   new TCanvas("cMuoHitP"   ,"mu",200,100); pMuoHitP->Draw();
 555   new TCanvas("cPioHitP"   ,"pi",200,100); pPioHitP->Draw();
 556   new TCanvas("cKaoHitP"   ,"K" ,200,100); pKaoHitP->Draw();
 557   new TCanvas("cProHitP"   ,"p" ,200,100); pProHitP->Draw();
 558   new TCanvas("cVertex"    ,"2d vertex" ,200,100); pVertex->Draw();
 559   new TCanvas("cRho"    ,"Rho of sec" ,200,100); pRho->Draw();
 560
 561   gBenchmark->Show("HitsPlots");
 562 }//HitQA()
 563
 564
 565
 566
 567