new way to calculate paritcle contribs
[u/mrichter/AliRoot.git] / RICH / menu.C
1 //globals for easy manual manipulations
2 AliRun *a;       
3 AliRunLoader *al;
4 AliLoader *rl;
5 AliRICH *r;
6 AliStack *s;
7
8 AliRICH * R()    {return r;}
9 void ph(Int_t event=0) {R()->PrintHits(event);}    //utility print hits for 'event' event
10 void ps(Int_t event=0) {R()->PrintSDigits(event);} //utility print sdigits
11 void pd(Int_t event=0) {R()->PrintDigits(event);}  //utility print digits
12 void pc(Int_t event=0) {R()->PrintClusters(event);}//utility print clusters
13 void pt(Int_t event=0) {R()->PrintTracks(event);}  //utility print tracks
14 //void st(Int_t event=0) {R()->PrintStack(event);} //utiltity print stack
15 //__________________________________________________________________________________________________
16 void pp(int tid)
17 {
18   if(!al) return;
19   al->LoadHeader();  al->LoadKinematics();
20   
21   if(tid<0||tid>=al->Stack()->GetNtrack())
22     cout<<"Valid tid number is 0-"<<al->Stack()->GetNtrack()-1<<" for this event.\n";
23   else
24     PrintParticleInfo(tid);
25   
26   al->UnloadKinematics();  al->UnloadHeader();
27 }
28 //__________________________________________________________________________________________________
29 void PrintParticleInfo(int tid)
30 {
31 // Prints particle info for a given TID
32   TParticle *p=al->Stack()->Particle(tid);
33   cout<<p->GetName()<<"("<<tid<<")";
34   if(p->GetMother(0)!=-1){cout<<" from "; PrintParticleInfo(p->GetMother(0));}
35   else                   {cout<<endl;} 
36 }    
37 //__________________________________________________________________________________________________
38 Int_t prim(Int_t tid)
39 {
40 // Provides mother TID for the given TID
41   R()->GetLoader()->GetRunLoader()->LoadHeader();  R()->GetLoader()->GetRunLoader()->LoadKinematics();
42   
43   if(tid<0||tid>=R()->GetLoader()->GetRunLoader()->Stack()->GetNtrack())
44     cout<<"Valid tid number is 0-"<<R()->GetLoader()->GetRunLoader()->Stack()->GetNtrack()-1<<" for this event.\n";
45   else
46     while(1){
47       TParticle *p=R()->GetLoader()->GetRunLoader()->GetAliRun()->Stack()->Particle(tid);
48       if(p->GetMother(0)==-1) break;
49       tid=p->GetMother(0);
50     }
51   
52   R()->GetLoader()->GetRunLoader()->UnloadKinematics();  R()->GetLoader()->GetRunLoader()->UnloadHeader();
53   return tid;
54 }
55
56 //__________________________________________________________________________________________________
57 void Show()
58 {  
59   Info("","\n\n\n");
60 //load all trees  
61   al->LoadHeader(); 
62     al->LoadKinematics();  
63       Bool_t isHits=!rl->LoadHits();  
64         Bool_t isSdigits=!rl->LoadSDigits();  
65           Bool_t isDigits=!rl->LoadDigits();//loaders
66             Bool_t isClusters=!rl->LoadRecPoints();
67   
68   
69   for(Int_t iEventN=0;iEventN<a->GetEventsPerRun();iEventN++){//events loop
70     Int_t iNparticles=al->Stack()->GetNtrack();
71     Int_t iNprims=al->Stack()->GetNprimary();
72     
73     Int_t iKPlusCounter=0,iKMinusCounter=0;    
74     for(Int_t iParticleN=0;iParticleN<iNparticles;iParticleN++){//stack loop
75       TParticle *pPart=al->Stack()->Particle(iParticleN);
76       switch(pPart->GetPdgCode()){
77         case kKPlus: iKPlusCounter++; break;
78         case kKMinus:iKMinusCounter++; break;
79       }
80     }//stack loop
81     
82     Info("Show-STA","Evt %i->   %i particles %i primaries  %i K+ %i K-",
83                      iEventN,   iNparticles,    iNprims,      iKPlusCounter,         iKMinusCounter);
84     
85     Int_t iHitsCounter=0;
86     Int_t iNentries=rl->TreeH()->GetEntries();
87     for(Int_t iEntryN=0;iEntryN<iNentries;iEntryN++){//TreeH loop
88       rl->TreeH()->GetEntry(iEntryN);//get current entry (prim)          
89       
90       for(Int_t iHitN=0;iHitN<r->Hits()->GetEntries();iHitN++){//hits loop
91         iHitsCounter++;
92         AliRICHhit *pHit = (AliRICHhit*)r->Hits()->At(iHitN);//get current hit
93         TParticle *pPart=al->Stack()->Particle(pHit->GetTrack());//get stack particle which produced the current hit
94         FillContribs(pPart->GetPdgCode(),pHit->C(),kFALSE);
95       }//hits loop
96       
97       if(iEntryN<7) Info("Show","Evt %i-> prim %4i has %4i hits from %s (,%7.2f,%7.2f)",
98                   iEventN,iEntryN, r->Hits()->GetEntries(), pPart->GetName(), pPart->Theta()*r2d,pPart->Phi()*r2d);
99     }//TreeH loop
100     Info("Show-HIT","Evt %i->   %i particles %i primaries  %i entries in TreeH %i hits",
101                      iEventN,   iNparticles,    iNprims,      iNentries,         iHitsCounter);
102     FillContribs(pPart->GetPdgCode(),pHit->C(),kTRUE);
103     
104     if(isSdigits){
105       rl->TreeS()->GetEntry(0);
106       Info("Show-SDI","Evt %i contains %5i sdigits",iEventN,r->SDigits()->GetEntries());
107     }
108     if(isDigits){
109       rl->TreeD()->GetEntry(0);
110       for(int i=1;i<=7;i++)
111         Info("Show-DIG","Evt %i chamber %i contains %5i digits",
112                                  iEventN,   i,           r->Digits(i)->GetEntries());
113     }else
114         Info("Show-DIG","There is no digits for this event");
115     if(isClusters){
116       rl->TreeR()->GetEntry(0);
117       for(int i=1;i<=7;i++)
118         Info("Show-CLU","Evt %i chamber %i contains %5i clusters",
119                                  iEventN,   i,           r->Clusters(i)->GetEntries());
120     }
121     cout<<endl;
122   }//events loop
123 //unload all trees    
124   rl->UnloadHits();  
125     if(isSdigits) rl->UnloadSDigits(); 
126       if(isDigits) rl->UnloadDigits(); 
127         if(isClusters) rl->UnloadRecPoints();
128           al->UnloadHeader();
129             al->UnloadKinematics();
130   
131   TVector counters=r->Counters();
132   
133   counters(9)=counters(0); for(Int_t i=1;i<=8;i++) counters(9)-=counters(i);
134   counters.Print();
135   cout<<endl;
136 }//void Show()
137 //__________________________________________________________________________________________________
138
139
140
141 Bool_t ReadAlice()
142 {
143   Info("ReadAlice","Tring to read ALICE from SIMULATED FILE.");
144   if(gAlice) delete gAlice;      
145   if(!(al=AliRunLoader::Open("galice.root","AlicE","update"))){
146     gSystem->Exec("rm -rf *.root *.dat");
147     Error("menu.C::ReadAlice","galice.root broken, removing all this garbage then init new one");
148     new AliRun("gAlice","Alice experiment system");
149     AliLog::SetModuleDebugLevel("RICH",1);
150     gAlice->Init("Config.C");
151     r=(AliRICH*)gAlice->GetDetector("RICH");
152     a=gAlice; //for manual convinience
153     return kFALSE;
154   }
155   al->LoadgAlice();//
156   if(!gAlice) Fatal("menu.C::ReadAlice","No gAlice in file");
157   a=al->GetAliRun();//provides pointer to AliRun object
158   a->SetDebug(0);    
159 //RICH      
160   if(!(r=(AliRICH*)gAlice->GetDetector("RICH"))) Warning("RICH/menu.C::ReadAlice","No RICH in file");
161   if(!(rl=al->GetLoader("RICHLoader")))          Warning("RICH/menu.C::ReadAlice","No RICH loader in file");        
162         
163   Info("ReadAlice","Run contains %i event(s)",gAlice->GetEventsPerRun());      
164   return kTRUE;
165 }
166 //__________________________________________________________________________________________________
167 void TestMenu()
168 {
169   TControlBar *pMenu = new TControlBar("vertical","RICH test");
170   pMenu->AddButton("Test segmentation"  ,"rp->TestSeg()"  ,"Test AliRICHParam segmentation methods");
171   pMenu->AddButton("Test response"      ,"rp->TestResp()" ,"Test AliRICHParam response methods");
172   pMenu->AddButton("Test transformation","rp->TestTrans()","Test AliRICHParam transformation methods");
173   pMenu->AddButton("Test opticals"      ,".x Opticals.h"  ,"Test optical properties");
174   pMenu->Show();  
175 }//TestMenu()
176 //__________________________________________________________________________________________________
177 void menu()
178
179   TControlBar *pMenu = new TControlBar("vertical","RICH main");
180        
181   if(ReadAlice()){//it's from file, show some info
182     pMenu->AddButton("Show",            "Show()",             "Shows the structure of events in files");
183     pMenu->AddButton("Display Fast",    "AliRICHDisplFast *d = new AliRICHDisplFast(); d->Exec();",        "Display Fast");
184     pMenu->AddButton("Control Plots",   "R()->ControlPlots()","Create some control histograms");
185     
186   }else{//it's aliroot, simulate
187     pMenu->AddButton("Debug ON",     "DebugON();",   "Switch debug on-off");   
188     pMenu->AddButton("Debug OFF",    "DebugOFF();",   "Switch debug on-off");   
189     pMenu->AddButton("Run",         "a->Run(1)",       "Process!");
190     pMenu->AddButton("Geo GUI",     "GeomGui()",       "Shows geometry"); 
191     pMenu->AddButton("Read RAW",    "ReadRaw()",       "Read a list of digits from test beam file"); 
192   }
193   pMenu->AddButton("Test submenu",    "TestMenu()",            "Shows test submenu");
194   pMenu->AddButton("Browser",         "new TBrowser;",         "Start ROOT TBrowser");
195   pMenu->AddButton("Quit",            ".q",                    "Close session");
196   pMenu->Show();
197 }//menu()
198 //__________________________________________________________________________________________________
199 void DebugOFF(){  Info("DebugOFF","");  AliLog::SetGlobalDebugLevel(0);}
200 void DebugON() {  Info("DebugON","");   AliLog::SetGlobalDebugLevel(AliLog::kDebug);}
201 //__________________________________________________________________________________________________
202 void ReadRaw()
203 {
204 // First interation for TB raw data reader  
205   char *sRawFileName="beam/run1822_4sigma.dat";
206   FILE *pRawFile=fopen(sRawFileName,"r");
207   if(!pRawFile){Error("ReadRaw","Something wrong with raw data file %s",sRawFileName); return;}
208   
209   Int_t iNpads=0,q=0,x=0,y=0;
210   
211   while(fscanf(pRawFile,"%i",&iNpads)){
212     Info("ReadRaw","%i pads fired",iNpads);
213     for(Int_t i=0;i<iNpads;i++)
214       fscanf(pRawFile,"%i %i %i",&q,&x,&y);
215   }
216   fclose(pRawFile);
217 }
218 //__________________________________________________________________________________________________
219 void FillContribs(Int_t part,Int_t C, Bool_t print)
220 {
221   static Int_t contrib[10][7][2];
222   static Bool_t kEnter=kFALSE;
223
224   C--; // chamber numbering from 1 to 7
225   if(!kEnter) {
226     for(Int_t i=0;i<10;i++) {
227       for(Int_t j=0;j<7;j++) {
228         for(Int_t k=0;k<2;k++) contrib[i][j][k]=0;
229       }
230     }
231     kEnter=kTRUE;
232   }
233
234   if(print) {
235     for(Int_t k=0;k<2;k++) {cout << "----------------Positives  to RICH ---------------" << endl;
236                             cout << " chamber    1    2    3     4     5     6    7    " << endl;
237                             cout << " -------------------------------------------------" << endl;
238       for(Int_t i=0;i<4;i++) {
239         if(i==0) cout << " e";
240         if(i==1) cout << "pi";
241         if(i==2) cout << " K";
242         if(i==3) cout << " p";
243         if(k==0) cout << "+         ";
244         if(k==1) cout << "-         ";
245         for(Int_t j=0;j<7;j++) {
246           cout << contrib[i][j][k] << "    ";
247         }
248           cout << endl;
249       }
250     }
251   }
252
253   // +ves
254   if(part==kPositron)    contrib[0][C][0]++;
255   if(part==kPiPlus)      contrib[1][C][0]++;
256   if(part==kKPlus)       contrib[2][C][0]++;
257   if(part==kProton)      contrib[3][C][0]++;
258
259   // -ves
260   if(part==kElectron)    contrib[0][C][1]++;
261   if(part==kPiMinus)     contrib[1][C][1]++;
262   if(part==kKMinus)      contrib[2][C][1]++;
263   if(part==kProtonBar)   contrib[3][C][1]++;
264 }
265 //__________________________________________________________________________________________________
266 void ParticleContribs()
267 {
268   Bool_t isHits    =!rl->LoadHits();
269   if(!isHits){Error("Exec","No hits. No contribs to calculate.");return;}
270   
271   TH1F *pDistH1 = new TH1F("distH1","length of charge track in amp gap;cm",100,0.,5.);
272
273   al->LoadHeader();  al->LoadKinematics();
274
275
276
277   for(Int_t iEventN=0;iEventN<gAlice->GetEventsPerRun();iEventN++){//events Loop
278     al->GetEvent(iEventN);
279     
280     AliStack *pStack=al->Stack();
281     Info("","Event %i-> Totaly %i primaries %i",iEventN,pStack->GetNtrack(),pStack->GetNprimary());
282     
283     
284     for(Int_t i=0;i<rl->TreeH()->GetEntries();i++){//prims loop
285       rl->TreeH()->GetEntry(i);
286       for(Int_t j=0;j<r->Hits()->GetEntries();j++){//hits loop for a given primary
287         AliRICHhit *pHit = (AliRICHhit*)r->Hits()->At(j);
288         TParticle *pParticle = al->Stack()->Particle(pHit->GetTrack());
289 //          if(pParticle->P()>1) FillContribs(pParticle->GetPdgCode(),pHit->C(),kFALSE);
290         if(pParticle->GetPDG()->Charge()!=0) pDistH1->Fill(pHit->Length());
291       }//hits loop
292     }//prims loop
293   }//events loop
294   FillContribs(0,0,kTRUE);
295   pDistH1->Draw();
296 }//ParticleContribs()
297 //__________________________________________________________________________________________________
298 void CheckPR()
299 {
300 //Pattern recognition wirh Stack particles
301   TFile *pFile = new TFile("$(HOME)/RPR.root","RECREATE","RICH Pattern Recognition");
302   TNtupleD *hn = new TNtupleD("hn","ntuple","Pmod:Charge:TrackTheta:TrackPhi:TrackX:TrackY:MinX:MinY:ChargeMIP:ThetaCerenkov:NPhotons:MipIndex");
303   printf("\n\n");
304     printf("Pattern Recognition done for event %5i",0);
305   for(Int_t iEvtN=0;iEvtN<R()->GetLoader()->GetRunLoader()->GetNumberOfEvents();iEvtN++) {
306     R()->GetLoader()->GetRunLoader()->GetEvent(iEvtN);
307     AliRICHTracker *tr = new AliRICHTracker();
308     tr->RecWithStack(hn);
309 //    Info("CheckPR","Pattern Recognition done for event %i \b",iEvtN);
310     printf("\b\b\b\b\b%5i",iEvtN+1);
311   }
312   printf("\n\n");
313   pFile->Write();pFile->Close();
314 }
315
316 //__________________________________________________________________________________________________
317 void GeomGui()
318 {
319   if(gGeoManager){ 
320     gGeoManager->GetTopVolume()->Draw(); 
321     AliRICHParam::ShowAxis();
322   }else 
323     new G3GeometryGUI;
324 }  
325