RICH/RichMenu.C

   1 //globals for easy manual manipulations
   2 AliRun *a;
   3 AliRunLoader *al;
   4 AliLoader *rl;
   5 AliRICH *r;
   6 AliStack *s;
   7
   8 void ph(Int_t event=0)  {r->PrintHits(event);}    //utility print hits for 'event' event
   9 void ps(Int_t event=0)  {r->PrintSDigits(event);} //utility print sdigits
  10 void pd(Int_t event=0)  {r->PrintDigits(event);}  //utility print digits
  11 void pc(Int_t event=0)  {r->PrintClusters(event);}//utility print clusters
  12 void pt(Int_t event=0)  {r->PrintTracks(event);}  //utility print tracks
  13 Int_t nem(Int_t event=0)  {AliRICHDisplFast::Nparticles(kElectron  ,event,al);} //utility number of electrons
  14 Int_t nep(Int_t event=0)  {AliRICHDisplFast::Nparticles(kPositron  ,event,al);} //utility number of positrons
  15 Int_t nmup(Int_t event=0) {AliRICHDisplFast::Nparticles(kMuonPlus  ,event,al);} //utility number of positive muons
  16 Int_t nmum(Int_t event=0) {AliRICHDisplFast::Nparticles(kMuonMinus ,event,al);} //utility number of negative muons
  17 Int_t npi0(Int_t event=0) {AliRICHDisplFast::Nparticles(kPi0       ,event,al);} //utility number of electrons
  18 Int_t npip(Int_t event=0) {AliRICHDisplFast::Nparticles(kPiPlus    ,event,al);} //utility number of electrons
  19 Int_t npim(Int_t event=0) {AliRICHDisplFast::Nparticles(kPiMinus   ,event,al);} //utility number of electrons
  20 Int_t nk0(Int_t event=0)  {AliRICHDisplFast::Nparticles(kK0        ,event,al);} //utility number of electrons
  21 Int_t nkp(Int_t event=0)  {AliRICHDisplFast::Nparticles(kKPlus     ,event,al);} //utility number of electrons
  22 Int_t nkm(Int_t event=0)  {AliRICHDisplFast::Nparticles(kKMinus    ,event,al);} //utility number of electrons
  23 Int_t npp(Int_t event=0)  {AliRICHDisplFast::Nparticles(kProton    ,event,al);} //utility number of protons
  24 //__________________________________________________________________________________________________
  25 void pp(int tid)
  26 {
  27   if(!al) return;
  28   al->LoadHeader();  al->LoadKinematics();
  29
  30   if(tid<0||tid>=al->Stack()->GetNtrack())
  31     cout<<"Valid tid number is 0-"<<al->Stack()->GetNtrack()-1<<" for this event.\n";
  32   else
  33     PrintParticleInfo(tid);
  34
  35   al->UnloadKinematics();  al->UnloadHeader();
  36 }
  37 //__________________________________________________________________________________________________
  38 void PrintParticleInfo(int tid)
  39 {
  40 // Prints particle info for a given TID
  41   TParticle *p=al->Stack()->Particle(tid);
  42   cout<<p->GetName()<<"("<<tid<<")";
  43   if(p->GetMother(0)!=-1){cout<<" from "; PrintParticleInfo(p->GetMother(0));}
  44   else                   {cout<<endl;}
  45 }
  46 //__________________________________________________________________________________________________
  47 Int_t prim(Int_t tid)
  48 {
  49 // Provides mother TID for the given TID
  50   al->LoadHeader();  al->LoadKinematics();
  51
  52   if(tid<0||tid>=al->Stack()->GetNtrack())
  53     cout<<"Valid tid number is 0-"<<al->Stack()->GetNtrack()-1<<" for this event.\n";
  54   else
  55     while(1){
  56       TParticle *p=al->Stack()->Particle(tid);
  57       if(p->GetMother(0)==-1) break;
  58       tid=p->GetMother(0);
  59     }
  60
  61   al->UnloadKinematics();  al->UnloadHeader();
  62   return tid;
  63 }
  64
  65 //__________________________________________________________________________________________________
  66 void Show()
  67 {
  68   CreateHists();
  69   Info("","\n\n\n");
  70 //load all trees
  71   al->LoadHeader();
  72     al->LoadKinematics();
  73       Bool_t isHits=!rl->LoadHits();
  74         Bool_t isSdigits=!rl->LoadSDigits();
  75           Bool_t isDigits=!rl->LoadDigits();//loaders
  76             Bool_t isClusters=!rl->LoadRecPoints();
  77
  78   for(Int_t iEventN=0;iEventN<a->GetEventsPerRun();iEventN++){//events loop
  79
  80
  81     Info("Show-STA","Evt %i->   %i particles %i primaries  %i e- %i e+",
  82                      iEventN,   iNparticles,    iNprims,      iElectronCounter,      iPositronCounter);
  83     Info("Show-STA","Evt %i->   %i particles %i primaries  %i p+ %i p-",
  84                      iEventN,   iNparticles,    iNprims,      iPiPlusCounter,        iPiMinusCounter);
  85     Info("Show-STA","Evt %i->   %i particles %i primaries  %i K+ %i K-",
  86                      iEventN,   iNparticles,    iNprims,      iKPlusCounter,         iKMinusCounter);
  87     Info("Show-STA","Evt %i->   %i particles %i primaries  %i p  %i pbar",
  88                      iEventN,   iNparticles,    iNprims,      iProtonCounter,        iProtonBarCounter);
  89
  90     Int_t iHitsCounter=0;
  91     Info("Show-HIT","Evt %i->   %i particles %i primaries  %i entries in TreeH %i hits",
  92                      iEventN,   iNparticles,    iNprims,      iNentries,         iHitsCounter);
  93
  94     if(isSdigits){
  95       rl->TreeS()->GetEntry(0);
  96       Info("Show-SDI","Evt %i contains %5i sdigits",iEventN,r->SDigits()->GetEntries());
  97     }
  98     if(isDigits){
  99       rl->TreeD()->GetEntry(0);
 100       for(int i=1;i<=7;i++)
 101         Info("Show-DIG","Evt %i chamber %i contains %5i digits",
 102                                  iEventN,   i,           r->Digits(i)->GetEntries());
 103     }else
 104         Info("Show-DIG","There is no digits for this event");
 105     if(isClusters){
 106       rl->TreeR()->GetEntry(0);
 107       for(int i=1;i<=7;i++)
 108         Info("Show-CLU","Evt %i chamber %i contains %5i clusters",
 109                                  iEventN,   i,           r->Clusters(i)->GetEntries());
 110     }
 111     cout<<endl;
 112   }//events loop
 113 //unload all trees
 114   rl->UnloadHits();
 115     if(isSdigits) rl->UnloadSDigits();
 116       if(isDigits) rl->UnloadDigits();
 117         if(isClusters) rl->UnloadRecPoints();
 118           al->UnloadHeader();
 119             al->UnloadKinematics();
 120   ShowHists();
 121 }//void Show()
 122 //__________________________________________________________________________________________________
 123
 124
 125
 126 Bool_t ReadAlice()
 127 {
 128   Info("ReadAlice","Tring to read ALICE from SIMULATED FILE.");
 129   if(gAlice){
 130     delete gAlice->GetRunLoader();
 131     delete gAlice;
 132   }
 133   if(!(al=AliRunLoader::Open())){//if not possible to read from galice.root, then create the new session
 134     gSystem->Exec("rm -rf *.root *.dat");
 135     Error("menu.C::ReadAlice","galice.root broken, removing all this garbage then init new one");
 136     new AliRun("gAlice","Alice experiment system");
 137     AliLog::SetModuleDebugLevel("RICH",1);
 138     gAlice->Init("Config.C");
 139     r=(AliRICH*)gAlice->GetDetector("RICH");
 140     a=gAlice; //for manual convinience
 141     return kFALSE;
 142   }
 143   al->LoadgAlice();//before this gAlice is 0;
 144   if(!gAlice) Fatal("menu.C::ReadAlice","No gAlice in file");
 145   a=al->GetAliRun();//provides pointer to AliRun object
 146 //RICH
 147   if(!(r=(AliRICH*)gAlice->GetDetector("RICH"))) Warning("RICH/menu.C::ReadAlice","No RICH in file");
 148   if(!(rl=al->GetLoader("RICHLoader")))          Warning("RICH/menu.C::ReadAlice","No RICH loader in file");
 149
 150   Info("ReadAlice","Run contains %i event(s)",gAlice->GetEventsPerRun());
 151   return kTRUE;
 152 }
 153 //__________________________________________________________________________________________________
 154 void TestMenu()
 155 {
 156   TControlBar *pMenu = new TControlBar("vertical","RICH test");
 157   pMenu->AddButton("Test segmentation"  ,"rp->TestSeg()"  ,"Test AliRICHParam segmentation methods");
 158   pMenu->AddButton("Test response"      ,"rp->TestResp()" ,"Test AliRICHParam response methods");
 159   pMenu->AddButton("Test transformation","rp->TestTrans()","Test AliRICHParam transformation methods");
 160   pMenu->AddButton("Test opticals"      ,".x Opticals.h"  ,"Test optical properties");
 161   pMenu->Show();
 162 }//TestMenu()
 163 //__________________________________________________________________________________________________
 164 void RichMenu()
 165 {
 166   TControlBar *pMenu = new TControlBar("vertical","RICH main");
 167
 168   if(ReadAlice()){//it's from file, show some info
 169     pMenu->AddButton("Show",            "Show()",             "Shows the structure of events in files");
 170     pMenu->AddButton("Display Fast",    "AliRICHDisplFast *d = new AliRICHDisplFast(); d->Exec();",        "Display Fast");
 171     pMenu->AddButton("Control Plots",   "ControlPlots()",     "Create some control histograms");
 172     pMenu->AddButton("Recon with stack","r->CheckPR()",     "Create RSR.root with ntuple hn");
 173   }else{//it's aliroot, simulate
 174     pMenu->AddButton("Debug ON",     "DebugON();",   "Switch debug on-off");
 175     pMenu->AddButton("Debug OFF",    "DebugOFF();",   "Switch debug on-off");
 176     pMenu->AddButton("Run",         "a->Run(1)",       "Process!");
 177     pMenu->AddButton("Geo GUI",     "GeomGui()",       "Shows geometry");
 178     pMenu->AddButton("Read RAW",    "ReadRaw()",       "Read a list of digits from test beam file");
 179   }
 180   pMenu->AddButton("Test submenu",    "TestMenu()",            "Shows test submenu");
 181   pMenu->AddButton("Browser",         "new TBrowser;",         "Start ROOT TBrowser");
 182   pMenu->AddButton("Quit",            ".q",                    "Close session");
 183   pMenu->Show();
 184 }//menu()
 185 //__________________________________________________________________________________________________
 186 void DebugOFF(){  Info("DebugOFF","");  AliLog::SetGlobalDebugLevel(0);}
 187 void DebugON() {  Info("DebugON","");   AliLog::SetGlobalDebugLevel(AliLog::kDebug);}
 188 //__________________________________________________________________________________________________
 189 void GeomGui()
 190 {
 191   if(gGeoManager){
 192     gGeoManager->GetTopVolume()->Draw();
 193     AliRICHParam::ShowAxis();
 194   }else
 195     new G3GeometryGUI;
 196 }
 197 //__________________________________________________________________________________________________
 198 void ControlPlots()
 199 {
 200   r->ControlPlots();
 201   new TBrowser;
 202 }
 203 //__________________________________________________________________________________________________
 204 void ParticleComposition()
 205 {
 206   TH2F *pFluxH2    =new TH2F("flux","Charged flux for central Hijing event with Vertex<470 P>5MeV for electrons and positrons, P>1GeV for others",10,-5,5, 10,0,10); pFluxH2->SetStats(0);
 207   pFluxH2->GetXaxis()->SetBinLabel(1,Form("p^{-}>%dGeV"   ,cutPproton));
 208   pFluxH2->GetXaxis()->SetBinLabel(2,Form("K^{-}>%dGeV"   ,cutPkaonminus));
 209   pFluxH2->GetXaxis()->SetBinLabel(3,Form("#pi^{-}>%dGeV" , ));
 210   pFluxH2->GetXaxis()->SetBinLabel(4,Form("#mu^{-}>%dGeV" , ));
 211   pFluxH2->GetXaxis()->SetBinLabel(5,Form("e^{+}>%dGeV"   ,));
 212
 213   pFluxH2->GetXaxis()->SetBinLabel(6,Form("e^{-}>%dGeV"));
 214   pFluxH2->GetXaxis()->SetBinLabel(7,Form("#mu^{+}>%dGeV"));
 215   pFluxH2->GetXaxis()->SetBinLabel(8,Form("#pi^{+}>%dGeV"));
 216   pFluxH2->GetXaxis()->SetBinLabel(9,Form("K^{+}>%dGeV"));
 217   pFluxH2->GetXaxis()->SetBinLabel(10,Form("p^{+}>%dGeV"));
 218
 219   pFluxH2->GetYaxis()->SetBinLabel(1,"sum");
 220   pFluxH2->GetYaxis()->SetBinLabel(2,"ch1");
 221   pFluxH2->GetYaxis()->SetBinLabel(3,"ch2");
 222   pFluxH2->GetYaxis()->SetBinLabel(4,"ch3");
 223   pFluxH2->GetYaxis()->SetBinLabel(5,"ch4");
 224   pFluxH2->GetYaxis()->SetBinLabel(6,"ch5");
 225   pFluxH2->GetYaxis()->SetBinLabel(7,"ch6");
 226   pFluxH2->GetYaxis()->SetBinLabel(8,"ch7");
 227   pFluxH2->GetYaxis()->SetBinLabel(9,"prim");
 228   pFluxH2->GetYaxis()->SetBinLabel(10,"tot");
 229
 230   TH1F *pElecP=new TH1F("Pelec","Electrons made hit in RICH;p [GeV]",1000,-10,10);
 231   TH1F *pMuonP=new TH1F("Pmuon","Muons made hit in RICH;p [GeV]"    ,1000,-10,10);
 232   TH1F *pPionP=new TH1F("Ppion","Pions made hit in RICH;p [GeV]"    ,1000,-10,10);
 233   TH1F *pKaonP=new TH1F("Pkaon","Kaon made hit in RICH;p [GeV]"     ,1000,-10,10);
 234   TH1F *pProtP=new TH1F("Pprot","Protons made hit in RICH;p [GeV]"  ,1000,-10,10);
 235
 236   al->LoadHeader();
 237   al->LoadKinematics();
 238   Int_t iNparticles=al->Stack()->GetNtrack();
 239   Int_t iNprims=al->Stack()->GetNprimary();
 240
 241   for(Int_t iParticleN=0;iParticleN<iNparticles;iParticleN++){//stack loop
 242     TParticle *pPart=al->Stack()->Particle(iParticleN);
 243
 244     if(iParticleN%10000==0) Info("Show-STA"," %i particles read",iParticleN);
 245
 246     switch(pPart->GetPdgCode()){
 247       case kPositron:  pFluxH2->Fill(-1,9);  break;
 248       case kElectron:  pFluxH2->Fill( 0,9);  break;
 249
 250       case kMuonMinus: pFluxH2->Fill(-2,9);  break;
 251       case kMuonPlus:  pFluxH2->Fill( 1,9);  break;
 252
 253       case kPiMinus:   pFluxH2->Fill(-3,9);  break;
 254       case kPiPlus:    pFluxH2->Fill( 2,9);  break;
 255
 256       case kKMinus:    pFluxH2->Fill(-4,9);  break;
 257       case kKPlus:     pFluxH2->Fill( 3,9);  break;
 258
 259       case kProtonBar: pFluxH2->Fill(-5,9);  break;
 260       case kProton:    pFluxH2->Fill( 4,9);  break;
 261     }//switch
 262   }//stack loop
 263
 264
 265
 266
 267   rl->LoadHits();
 268
 269
 270   for(Int_t iEntryN=0;iEntryN < rl->TreeH()->GetEntries();iEntryN++){//TreeH loop
 271     rl->TreeH()->GetEntry(iEntryN);//get current entry (prim)
 272     for(Int_t iHitN=0;iHitN<r->Hits()->GetEntries();iHitN++){//hits loop
 273       AliRICHhit *pHit = (AliRICHhit*)r->Hits()->At(iHitN);//get current hit
 274       TParticle  *pPart=al->Stack()->Particle(pHit->GetTrack());//get stack particle which produced the current hit
 275
 276       if(TMath::Sqrt(pPart->Vx()*pPart->Vx()+pPart->Vy()*pPart->Vy()+pPart->Vz()*pPart->Vz()) > 470) continue; //cut on vertex position
 277       if(pPart->GetPdgCode()==kElectron && pPart->P()<0.005)                                         continue; //cut on electron momentum 5MeV
 278       if(pPart->GetPdgCode()==kPositron && pPart->P()<0.005)                                         continue; //cut on electron momentum 5MeV
 279       if(pPart->GetPdgCode()!=kElectron && pPart->P()<1)                                             continue; //cut on others momentum 1GeV
 280
 281       switch(pPart->GetPdgCode()){
 282         case kPositron : pElecP->Fill(-pPart->P()); pFluxH2->Fill(-1,pHit->C());break;
 283         case kElectron : pElecP->Fill( pPart->P()); pFluxH2->Fill( 0,pHit->C());break;
 284
 285         case kMuonPlus : pMuonP->Fill( pPart->P()); pFluxH2->Fill(-2,pHit->C());break;
 286         case kMuonMinus: pMuonP->Fill(-pPart->P()); pFluxH2->Fill( 1,pHit->C());break;
 287
 288         case kPiMinus  : pPionP->Fill(-pPart->P()); pFluxH2->Fill(-3,pHit->C());break;
 289         case kPiPlus   : pPionP->Fill( pPart->P()); pFluxH2->Fill( 2,pHit->C());break;
 290
 291         case kKPlus    : pKaonP->Fill( pPart->P()); pFluxH2->Fill(-4,pHit->C());break;
 292         case kKMinus   : pKaonP->Fill(-pPart->P()); pFluxH2->Fill( 3,pHit->C());break;
 293
 294         case kProton   : pProtP->Fill( pPart->P()); pFluxH2->Fill(-5,pHit->C());break;
 295         case kProtonBar: pProtP->Fill(-pPart->P()); pFluxH2->Fill( 4,pHit->C());break;
 296       }
 297     }//hits loop
 298   }//TreeH loop
 299
 300   rl->UnloadHits();
 301   al->UnloadHeader();
 302   al->UnloadKinematics();
 303
 304
 305   for(Int_t i=1;i<=pFluxH2->GetNbinsX();i++){
 306     Stat_t sum=0;
 307     for(Int_t j=2;j<=8;j++)    sum+=pFluxH2->GetBinContent(i,j);
 308     pFluxH2->SetBinContent(i,1,sum);
 309   }
 310
 311   TCanvas *pC1=new TCanvas("canvas1",Form("Event Nprims=%i",iNprims),1000,900);
 312   pFluxH2->Draw("text");  gPad->SetGrid();
 313
 314   new TCanvas("celec","",200,100); pElecP->Draw();
 315   new TCanvas("cmuon","",200,100); pMuonP->Draw();
 316   new TCanvas("cpion","",200,100); pPionP->Draw();
 317   new TCanvas("ckaon","",200,100); pKaonP->Draw();
 318   new TCanvas("cprot","",200,100); pProtP->Draw();
 319
 320 }
 321
 322
 323
 324
 325
 326
 327
 328
 329