]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HMPID/Hmenu.C
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Do not recreate existing containers
[u/mrichter/AliRoot.git] / HMPID / Hmenu.C
CommitLineData
594d861b 1AliRun *a; AliRunLoader *al; TGeoManager *g; //globals for easy manual manipulations
2AliHMPID *r; AliLoader *rl; AliHMPIDParam *rp;
3Bool_t isGeomType=kFALSE;
4
5//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
6void GetParam()
7{
8 isGeomType=!isGeomType;
9 if(g) delete g; if(rp) delete rp; //delete current TGeoManager and AliHMPIDParam
10 if(isGeomType) g=TGeoManager::Import("geometry.root");
11 else g=TGeoManager::Import("misaligned_geometry.root");
12 rp=AliHMPIDParam::Instance();
13}
14//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
15void Hmenu()
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("HMPID"); r=(AliHMPID*)a->GetDetector("HMPID"); //get HMPID object from galice.root
24
25 status+=Form(" with %i event(s)",al->GetNumberOfEvents()); status+=(r)? " with HMPID": " without HMPID";
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//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
45void 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//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
55void 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//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
69void 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//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
82void 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" ,"AliHMPIDParam::TestResp();","Test AliHMPIDParam response methods" );
88 pMenu->AddButton("Print map" ,"PrintMap();" ,"Test AliHMPIDParam transformation methods" );
89 pMenu->AddButton("Test Recon" ,"rec();" ,"Test AliHMPIDRecon" );
90 pMenu->Show();
91}//Test()
92//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
93
94
95void doff(){ Printf("DebugOFF"); AliLog::SetGlobalDebugLevel(0);}
96void don() { Printf("DebugON"); AliLog::SetGlobalDebugLevel(AliLog::kDebug);}
97
98void geo ( ) {gGeoManager->GetTopVolume()->Draw("ogl");}
99
100void du ( ) {r->Dump ( );} //utility display
101
102void hp (Int_t evt=0 ) {r->HitPrint (evt);} //print hits for requested event
103void hq ( ) {r->HitQA ( );} //hits QA plots for all events
104void sp (Int_t evt=0 ) {r->SdiPrint (evt);} //print sdigits for requested event
105void sq (Int_t evt=0 ) {r->SdiPrint (evt);} //print sdigits for requested event
106void dp (Int_t evt=0 ) {r->DigPrint (evt);} //print digits for requested event
107void dq ( ) {AliHMPIDReconstructor::DigQA (al );} //digits QA plots for all events
108void cp (Int_t evt=0 ) {r->CluPrint (evt); } //print clusters for requested event
109void cq ( ) {AliHMPIDReconstructor::CluQA (al );} //clusters QA plots for all events
110
111void ep ( ) {AliHMPIDTracker::EsdQA(1); }
112void eq ( ) {AliHMPIDTracker::EsdQA(); }
113void mp (Double_t probCut=0.7 ) {AliHMPIDTracker::MatrixPrint(probCut);}
114
115
116void t (Int_t evt=0 ) {AliHMPIDParam::Stack(evt);}
117void tid (Int_t tid,Int_t evt=0) {AliHMPIDParam::Stack(evt,tid);}
118//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
119void tst()
120{
121//create all lists
122 TClonesArray hit("AliHMPIDHit"),sdi("AliHMPIDDigit");
123 TObjArray dig,clu,cog; for(Int_t i=0;i<7;i++) {dig.AddAt(new TClonesArray("AliHMPIDDigit"),i);
124 clu.AddAt(new TClonesArray("AliHMPIDCluster"),i);
125 cog.AddAt(new TClonesArray("AliHMPIDCluster"),i);}
126//simulate track
127 TLorentzVector mom; mom.SetPtEtaPhiM(3,0,30*TMath::DegToRad(),0.938);
128 TLorentzVector vtx; vtx.SetXYZT(0,0,0,0);
129
130 AliESD *pESD=new AliESD; pESD->SetMagneticField(0.2);
131 pESD->AddTrack(new AliESDtrack(new TParticle(kProton,0,-1,-1,-1,-1,mom,vtx)));
132 pESD->Print();
133
134 AliTracker::SetFieldMap(new AliMagF,1);
135 AliHMPIDTracker *pTracker=new AliHMPIDTracker;
136
137
138 return;
139 Double_t th=8*TMath::DegToRad(); //gRandom->Rndm()*TMath::PiOver4();
140 Double_t ph=gRandom->Rndm()*TMath::TwoPi();
141 Double_t radx=gRandom->Rndm()*AliHMPIDDigit::SizeAllX();
142 Double_t rady=gRandom->Rndm()*AliHMPIDDigit::SizeAllY();
143
144 Int_t iHitCnt=0;
145
146
147
148
149
150
151 AliHMPIDv1::Hit2Sdi(&hitLst,&sdiLst);
152 AliHMPIDDigitizer::Sdi2Dig(&sdiLst,&digLst);
153 AliHMPIDReconstructor::Dig2Clu(&digLst,&cluLst); AliHMPIDReconstructor::Dig2Clu(&digLst,&cogLst,0);
154
155
156 Int_t iDigN=0,iCluN=0,iCogN=0; for(Int_t i=0;i<7;i++){ iDigN+=((TClonesArray*)digLst.At(i))->GetEntries();
157 iCluN+=((TClonesArray*)cluLst.At(i))->GetEntries();
158 iCogN+=((TClonesArray*)cogLst.At(i))->GetEntries(); }
159
160 Printf("SDI------SDI---------SDI--------SDI------SDI------SDI #%i",sdiLst.GetEntries());sdiLst.Print();Printf("");
161 Printf("DIG------DIG---------DIG--------DIG------DIG------DIG #%i",iDigN );digLst.Print();Printf("");
162 Printf("HIT------HIT---------HIT--------HIT------HIT------HIT #%i",hitLst.GetEntries());hitLst.Print();Printf("");
163 Printf("CLU------CLU---------CLU--------CLU------CLU------CLU #%i",iCluN );cluLst.Print();Printf("");
164 Printf("COG------COG---------COG--------COG------COG------COG #%i",iCogN );cogLst.Print();Printf("");
165}
166
167
168
169void PrintMap()
170{
171
172 Double_t r2d=TMath::RadToDeg();
173
174 Double_t x=AliHMPIDDigit::SizeAllX(),y=AliHMPIDDigit::SizeAllY();
175
176 Printf("\n\n\n");
177
178 for(int ch=6;ch>=0;ch--){
179 AliHMPIDDigit dL(ch,0,1,0,0),dR(ch,0,1,67,0);
180 TVector3 lt=rp->Lors2Mars(ch,0,y); TVector3 rt=rp->Lors2Mars(ch,x,y);
181 TVector3 ce=rp->Lors2Mars(ch,x/2,y/2);
182 TVector3 lb=rp->Lors2Mars(ch,0,0); TVector3 rb=rp->Lors2Mars(ch,x,0);
183
184 Printf(" ____________________________");
185 Printf("|%6.2fcm %6.2fcm|" ,lt.Mag() , rt.Mag() );
186 Printf("|%6.2fde %6.2fde|" ,lt.Theta()*r2d , rt.Theta()*r2d );
187 Printf("|%6.2fde %6.2fde|" ,lt.Phi()*r2d , rt.Phi()*r2d );
188 Printf("| |" );
189 Printf("|DDL %2i %7.2fcm DDL %2i|" ,dL.DdlIdx() , ce.Mag() , dR.DdlIdx() );
190 Printf("| 0x%x %7.2fdeg 0x%x|" ,dL.DdlId() , ce.Theta()*r2d , dR.DdlId() );
191 Printf("| %7.2fdeg |" , ce.Phi()*r2d );
192 Printf("| |");
193 Printf("|%6.2fcm %6.2fcm|" ,lb.Mag() , rb.Mag() );
194 Printf("|%6.2fde %6.2fde|" ,lb.Theta()*r2d , rb.Theta()*r2d );
195 Printf("|%6.2fde Ch%i %6.2fde|" ,lb.Phi()*r2d , ch , rb.Phi()*r2d );
196 Printf(" ----------------------------");
197 }
198
199 Double_t m[3];
200 for(int i=0;i<1000;i++){
201 Float_t xout=0,xin=gRandom->Rndm()*130.60;
202 Float_t yout=0,yin=gRandom->Rndm()*126.16;
203 Int_t c=gRandom->Rndm()*6;
204 rp->Lors2Mars(c,xin,yin,m);
205 rp->Mars2Lors(c,m,xout,yout);
206 if( (xin-xout) != 0) Printf("Problem in X");
207 if( (yin-yout) != 0) Printf("Problem in Y");
208 }
209
210 Int_t ddl,r,d,a,ch,raw,pc,px,py; AliHMPIDDigit dd;
211
212 ddl=0;raw=0x2214000;r= 8;d=8;a=20;
213 ddl=1;raw=0x2214000;r= 8;d=8;a=20;
214
215
216 ddl=2;raw=0x08d6000;r= 2;d=3;a=22;
217 ddl=3;raw=0x08d6000;r= 2;d=3;a=22;
218
219
220 ddl=6;raw=0x592e000;r=22;d=4;a=46;ch=3;pc=4;px=55;py=5;dd.Raw(ddl,raw);
221 Printf("(ch=%i,pc=%i,x=%2i,y=%2i) ddl=%i raw=0x%h (r=%2i,d=%2i,a=%2i)",
222 ch, pc, px, py, ddl, raw, r, d, a); dd.Print();
223 ddl=7;raw=0x592e000;r=22;d=4;a=46;ch=3;pc=1;
224
225
226}//PrintMap()
227
228
229void rec()
230{
231
232
233 Double_t ckovMax=0.75,ckovSim;
234 Int_t nSim=0;
235 while(nSim<3){
236 ckovSim=gRandom->Rndm()*ckovMax;//0.6468;
237 nSim=20*TMath::Power(TMath::Sin(ckovSim)/TMath::Sin(ckovMax),2); //scale number of photons
238 }
239
240
241 TClonesArray *pCluLst=new TClonesArray("AliHMPIDCluster");
242 TPolyMarker *pMipMap=new TPolyMarker(); pMipMap->SetMarkerStyle(8); pMipMap->SetMarkerColor(kRed);
243 TPolyMarker *pPhoMap=new TPolyMarker(); pPhoMap->SetMarkerStyle(4); pPhoMap->SetMarkerColor(kRed);
244 TPolyMarker *pBkgMap=new TPolyMarker(); pBkgMap->SetMarkerStyle(25); pBkgMap->SetMarkerColor(kRed);
245 TPolyLine *pRing =new TPolyLine; pRing->SetLineColor(kGreen);
246 TPolyLine *pOld =new TPolyLine; pOld->SetLineColor(kBlue);
247
248 Int_t iCluCnt=0; pMipMap->SetPoint(iCluCnt,x,y); new((*pCluLst)[iCluCnt++]) AliHMPIDCluster(1,x,y,200); //add mip cluster
249
250// for(int j=0;j<30;j++){ //add bkg photons
251// Float_t bkgX=gRandom->Rndm()*AliHMPIDDigit::SizeAllX();
252// Float_t bkgY=gRandom->Rndm()*AliHMPIDDigit::SizeAllY();
253// pBkgMap->SetPoint(iCluCnt,bkgX,bkgY); new((*pCluLst)[iCluCnt++]) AliHMPIDCluster(1,bkgX,bkgY,35);
254// }
255
256
257
258
259 AliHMPIDRecon rec; rec.SetTrack(th,ph,x,y);
260
261 TVector2 pos;
262 for(int i=0;i<nSim;i++){
263 rec.TracePhot(ckovSim,gRandom->Rndm()*2*TMath::Pi(),pos); //add photons
264 if(AliHMPIDDigit::IsInDead(pos.X(),pos.Y())) continue;
265 pPhoMap->SetPoint(iCluCnt,pos.X(),pos.Y()); new((*pCluLst)[iCluCnt++]) AliHMPIDCluster(1,pos.X(),pos.Y(),35);
266 }
267
268
269 Int_t nRec=0,nOld=0;
270 Double_t ckovRec=rec.CkovAngle(pCluLst,nRec); Double_t err=TMath::Sqrt(rec.CkovSigma2());
271 Double_t ckovOld=old.CkovAngle(pCluLst,nOld);
272
273 Printf("---------------- Now reconstructed --------------------");
274
275
276 for(int j=0;j<100;j++){rec.TracePhot(ckovRec,j*0.0628,pos); pRing->SetPoint(j,pos.X(),pos.Y());}
277 for(int j=0;j<100;j++){rec.TracePhot(ckovOld,j*0.0628,pos); pOld->SetPoint(j,pos.X(),pos.Y());}
278
279 new TCanvas; AliHMPIDDigit::DrawPc(); pMipMap->Draw(); pPhoMap->Draw(); pBkgMap->Draw(); pRing->Draw(); pOld->Draw();
280
281 TLatex txt; txt.SetTextSize(0.03);
282 txt.DrawLatex(65,127,Form("#theta=%.4f#pm%.5f with %2i #check{C}" ,ckovSim, 0.,nSim ));
283 txt.DrawLatex(65,122,Form("#theta=%.4f#pm%.5f with %2i #check{C} Old=%.4f with %i #check{C}" ,ckovRec,err,nRec,ckovOld,nOld));
284 txt.DrawLatex(0 ,127,Form("#theta=%.2f#circ #phi=%.2f#circ @(%.2f,%.2f) ",th*TMath::RadToDeg(),ph*TMath::RadToDeg(),x,y));
285
286// for(int i=0;i<35;i++){
287// Double_t ckov=0.1+i*0.02;
288// Printf("Ckov=%.2f Old=%.3f New=%.3f",ckov,old.FindRingArea(ckov),rec.FindRingArea(ckov));
289// }
290
291
292}//rec()
293
294
295void HitQA(Double_t cut,Double_t cutele,Double_t cutR)
296{
297// Provides a set of control plots intended primarily for charged particle flux analisys
298// Arguments: cut (GeV) - cut on momentum of any charged particles but electrons,
299// cetele (GeV) - the same for electrons-positrons
300// cutR (cm) - cut on production vertex radius (cylindrical system)
301 gBenchmark->Start("HitsAna");
302
303 Double_t cutPantiproton =cut;
304 Double_t cutPkaonminus =cut;
305 Double_t cutPpionminus =cut;
306 Double_t cutPmuonminus =cut;
307 Double_t cutPpositron =cutele;
308
309 Double_t cutPelectron =cutele;
310 Double_t cutPmuonplus =cut;
311 Double_t cutPpionplus =cut;
312 Double_t cutPkaonplus =cut;
313 Double_t cutPproton =cut;
314
315
316 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
317 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
318 TH1F *pEleAllP =new TH1F("EleAllP" , "e^{+} e^{-} all;p[GeV]" ,1000,-1 ,1 );
319 TH1F *pEleHitP =new TH1F("EleHitP" ,Form("e^{+} e^{-} hit %s;p[GeV]" ,GetName()) ,1000,-1 ,1 );
320 TH1F *pMuoHitP =new TH1F("MuoHitP" ,Form("#mu^{-} #mu^{+} hit %s;p[GeV]" ,GetName()) ,1000,-4 ,4 );
321 TH1F *pPioHitP =new TH1F("PioHitP" ,Form("#pi^{-} #pi^{+} hit %s;p[GeV]" ,GetName()) ,1000,-4 ,4 );
322 TH1F *pKaoHitP =new TH1F("KaoHitP" ,Form("K^{-} K^{+} hit %s;p[GeV]" ,GetName()) ,1000,-4 ,4 );
323 TH1F *pProHitP =new TH1F("ProHitP" ,Form("p^{-} p^{+} hit %s;p[GeV]" ,GetName()) ,1000,-4 ,4 );
324 TH2F *pFlux =new TH2F("flux" ,Form("%s flux with Rvertex<%.1fcm" ,GetName(),cutR),10 ,-5 ,5 , 10,0 ,10); //special text hist
325 TH2F *pVertex =new TH2F("vertex" ,Form("%s 2D vertex of HMPID hit;x;y" ,GetName()) ,120 ,0 ,600 ,120,0 ,600); //special text hist
326 TH1F *pRho =new TH1F("rho" ,Form("%s r of HMPID hit" ,GetName()) ,600 ,0 ,600); //special text hist
327 pFlux->SetStats(0);
328 pFlux->GetXaxis()->SetBinLabel(1 ,Form("p^{-}>%.3fGeV/c" ,cutPantiproton));
329 pFlux->GetXaxis()->SetBinLabel(2 ,Form("K^{-}>%.3fGeV/c" ,cutPkaonminus ));
330 pFlux->GetXaxis()->SetBinLabel(3 ,Form("#pi^{-}>%.3fGeV/c" ,cutPpionminus ));
331 pFlux->GetXaxis()->SetBinLabel(4 ,Form("#mu^{-}>%.3fGeV/c" ,cutPmuonminus ));
332 pFlux->GetXaxis()->SetBinLabel(5 ,Form("e^{+}>%.3fGeV/c" ,cutPpositron ));
333
334 pFlux->GetXaxis()->SetBinLabel(6 ,Form("e^{-}>%.3fGeV/c" ,cutPelectron ));
335 pFlux->GetXaxis()->SetBinLabel(7 ,Form("#mu^{+}>%.3fGeV/c" ,cutPmuonplus ));
336 pFlux->GetXaxis()->SetBinLabel(8 ,Form("#pi^{+}>%.3fGeV/c" ,cutPpionplus ));
337 pFlux->GetXaxis()->SetBinLabel(9 ,Form("K^{+}>%.3fGeV/c" ,cutPkaonplus ));
338 pFlux->GetXaxis()->SetBinLabel(10,Form("p^{+}>%.3fGeV/c" ,cutPproton ));
339
340 pFlux->GetYaxis()->SetBinLabel(1,"sum");
341 pFlux->GetYaxis()->SetBinLabel(2,"ch1");
342 pFlux->GetYaxis()->SetBinLabel(3,"ch2");
343 pFlux->GetYaxis()->SetBinLabel(4,"ch3");
344 pFlux->GetYaxis()->SetBinLabel(5,"ch4");
345 pFlux->GetYaxis()->SetBinLabel(6,"ch5");
346 pFlux->GetYaxis()->SetBinLabel(7,"ch6");
347 pFlux->GetYaxis()->SetBinLabel(8,"ch7");
348 pFlux->GetYaxis()->SetBinLabel(9,"prim");
349 pFlux->GetYaxis()->SetBinLabel(10,"tot");
350
351//end of hists definition
352
353 Int_t iNevents=fLoader->GetRunLoader()->GetAliRun()->GetEventsPerRun(),iCntPrimParts=0,iCntTotParts=0;
354//load all needed trees
355 fLoader->LoadHits();
356 fLoader->GetRunLoader()->LoadHeader();
357 fLoader->GetRunLoader()->LoadKinematics();
358
359 for(Int_t iEvtN=0;iEvtN < iNevents;iEvtN++){//events loop
360 fLoader->GetRunLoader()->GetEvent(iEvtN);
361 AliInfo(Form(" %i event processes",fLoader->GetRunLoader()->GetEventNumber()));
362 AliStack *pStack= fLoader->GetRunLoader()->Stack();
363
364 for(Int_t iParticleN=0;iParticleN<pStack->GetNtrack();iParticleN++){//stack loop
365 TParticle *pPart=pStack->Particle(iParticleN);
366
367 if(iParticleN%10000==0) AliInfo(Form(" %i particles read",iParticleN));
368
369 switch(pPart->GetPdgCode()){
370 case kProtonBar: pFlux->Fill(-4.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-4.5,8); break;
371 case kKMinus: pFlux->Fill(-3.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-3.5,8); break;
372 case kPiMinus: pFlux->Fill(-2.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-2.5,8); break;
373 case kMuonMinus: pFlux->Fill(-1.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-1.5,8); break;
374 case kPositron: pFlux->Fill(-0.5,9); if(pPart->Rho()<0.01) pFlux->Fill(-0.5,8); pEleAllP->Fill(-pPart->P()); break;
375
376 case kElectron: pFlux->Fill( 0.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 0.5,8); pEleAllP->Fill( pPart->P()); break;
377 case kMuonPlus: pFlux->Fill( 1.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 1.5,8); break;
378 case kPiPlus: pFlux->Fill( 2.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 2.5,8); break;
379 case kKPlus: pFlux->Fill( 3.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 3.5,8); break;
380 case kProton: pFlux->Fill( 4.5,9); if(pPart->Rho()<0.01) pFlux->Fill( 4.5,8); break;
381 }//switch
382 }//stack loop
383//now hits analiser
384 for(Int_t iEntryN=0;iEntryN < fLoader->TreeH()->GetEntries();iEntryN++){//TreeH loop
385 fLoader->TreeH()->GetEntry(iEntryN); //get current entry (prim)
386 for(Int_t iHitN=0;iHitN < Hits()->GetEntries();iHitN++){//hits loop
387 AliHMPIDHit *pHit = (AliHMPIDHit*)Hits()->At(iHitN); //get current hit
388 TParticle *pPart=pStack->Particle(pHit->GetTrack()); //get stack particle which produced the current hit
389
390 if(pPart->GetPDG()->Charge()!=0&&pPart->Rho()>0.1) pVertex->Fill(pPart->Vx(),pPart->Vy()); //safe margin for sec.
391 if(pPart->GetPDG()->Charge()!=0) pRho->Fill(pPart->Rho()); //safe margin for sec.
392 if(pPart->R()>cutR) continue; //cut on production radius (cylindrical system)
393
394 switch(pPart->GetPdgCode()){
395 case kProtonBar: if(pPart->P()>cutPantiproton) {pProHitP->Fill(-pPart->P()); pFlux->Fill(-4.5,pHit->Ch());}break;
396 case kKMinus : if(pPart->P()>cutPkaonminus) {pKaoHitP->Fill(-pPart->P()); pFlux->Fill(-3.5,pHit->Ch());}break;
397 case kPiMinus : if(pPart->P()>cutPpionminus) {pPioHitP->Fill(-pPart->P()); pFlux->Fill(-2.5,pHit->Ch());}break;
398 case kMuonMinus: if(pPart->P()>cutPmuonminus) {pMuoHitP->Fill(-pPart->P()); pFlux->Fill(-1.5,pHit->Ch());}break;
399 case kPositron : if(pPart->P()>cutPpositron) {pEleHitP->Fill(-pPart->P()); pFlux->Fill(-0.5,pHit->Ch());
400 pEleHitRP->Fill(-pPart->P(),pPart->R()); pEleHitRZ->Fill(pPart->Vz(),pPart->R()); }break;
401
402 case kElectron : if(pPart->P()>cutPelectron) {pEleHitP->Fill( pPart->P()); pFlux->Fill( 0.5,pHit->Ch());
403 pEleHitRP->Fill( pPart->P(),pPart->R()); pEleHitRZ->Fill(pPart->Vz(),pPart->R()); }break;
404 case kMuonPlus : if(pPart->P()>cutPmuonplus) {pMuoHitP->Fill( pPart->P()); pFlux->Fill( 1.5,pHit->Ch());}break;
405 case kPiPlus : if(pPart->P()>cutPpionplus) {pPioHitP->Fill( pPart->P()); pFlux->Fill( 2.5,pHit->Ch());}break;
406 case kKPlus : if(pPart->P()>cutPkaonplus) {pKaoHitP->Fill( pPart->P()); pFlux->Fill( 3.5,pHit->Ch());}break;
407 case kProton : if(pPart->P()>cutPproton) {pProHitP->Fill( pPart->P()); pFlux->Fill( 4.5,pHit->Ch());}break;
408 }
409 }//hits loop
410 }//TreeH loop
411 iCntPrimParts +=pStack->GetNprimary();
412 iCntTotParts +=pStack->GetNtrack();
413 }//events loop
414//unload all loaded staff
415 fLoader->UnloadHits();
416 fLoader->GetRunLoader()->UnloadHeader();
417 fLoader->GetRunLoader()->UnloadKinematics();
418//Calculater some sums
419 Stat_t sum=0;
420//sum row, sum over rows
421 for(Int_t i=1;i<=pFlux->GetNbinsX();i++){
422 sum=0; for(Int_t j=2;j<=8;j++) sum+=pFlux->GetBinContent(i,j);
423 pFlux->SetBinContent(i,1,sum);
424 }
425
426//display everything
427 new TCanvas("canvas1",Form("Events %i Nprims=%i Nparticles=%i",iNevents,iCntPrimParts,iCntTotParts),1000,900); pFlux->Draw("text"); gPad->SetGrid();
428//total prims and particles
429 TLatex latex; latex.SetTextSize(0.02);
430 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));
431 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));
432 for(Int_t iCh=0;iCh<7;iCh++) {
433 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));
434 }
435 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));
436
437 new TCanvas("cEleAllP" ,"e" ,200,100); pEleAllP->Draw();
438 new TCanvas("cEleHitRP" ,"e" ,200,100); pEleHitRP->Draw();
439 new TCanvas("cEleHitRZ" ,"e" ,200,100); pEleHitRZ->Draw();
440 new TCanvas("cEleHitP" ,"e" ,200,100); pEleHitP->Draw();
441 new TCanvas("cMuoHitP" ,"mu",200,100); pMuoHitP->Draw();
442 new TCanvas("cPioHitP" ,"pi",200,100); pPioHitP->Draw();
443 new TCanvas("cKaoHitP" ,"K" ,200,100); pKaoHitP->Draw();
444 new TCanvas("cProHitP" ,"p" ,200,100); pProHitP->Draw();
445 new TCanvas("cVertex" ,"2d vertex" ,200,100); pVertex->Draw();
446 new TCanvas("cRho" ,"Rho of sec" ,200,100); pRho->Draw();
447
448 gBenchmark->Show("HitsPlots");
449}//HitQA()
450
451
452void RecWithStack()
453{
454 al->LoadHeader();al->LoadKinematics();
455 AliStack *pStk=al->Stack();
456
457 AliESD *pEsd=new AliESD;
458 for(Int_t iTrk=0;iTrk<pStk->GetNtrack();iTrk++){//stack loop
459 TParticle *pPart=pStk->Particle(iTrk);
460 if(pPart->GetPDG()->Charge()==0) continue; //neutral particles are not reconstructed
461 pEsd->AddTrack(new AliESDtrack(pPart));
462 }//stack loop
463
464 pEsd->Print();
465 AliTracker::SetFieldMap(new AliMagF,1);
466 AliHMPIDTracker t;
467 rl->LoadRecPoints();
468 t.LoadClusters(rl->TreeR());
469 t.PropagateBack(pEsd);
470 rl->UnloadRecPoints();
471}
472
473
474void AliHMPIDReconstructor::CluQA(AliRunLoader *pAL)
475{
476// Quality assesment plots for clusters.
477// This methode takes list of digits and form list of clusters again in order to
478// calculate cluster shape and cluster particle mixture
479 AliLoader *pRL=pAL->GetDetectorLoader("HMPID"); AliHMPID *pRich=(AliHMPID*)pAL->GetAliRun()->GetDetector("HMPID");//get pointers for HMPID and HMPID loader
480 Int_t iNevt=pAL->GetNumberOfEvents(); if(iNevt==0) {AliInfoClass("No events");return;}
481 if(pRL->LoadDigits()) {AliInfoClass("No digits file");return;}
482 pAL->LoadHeader();
483 pAL->LoadKinematics();
484// AliStack *pStack=pAL->Stack();
485 TH1::AddDirectory(kFALSE);
486
487
488 TH1F* pQ=new TH1F("RiAllQ" ,"Charge All" ,4000 ,0 ,4000);// Q hists
489 TH1F* pCerQ=new TH1F("RiCerQ" ,"Charge Ckov" ,4000 ,0 ,4000);
490 TH1F* pMipQ=new TH1F("RiMipQ" ,"Charge MIP" ,4000 ,0 ,4000);
491
492 TH1F* pS=new TH1F("RichCluSize" ,"Cluster size;size" ,100 ,0 ,100 );// size hists
493 TH1F* pCerS=new TH1F("RichCluCerSize" ,"Ckov size;size" ,100 ,0 ,100 );
494 TH1F* pMipS=new TH1F("RichCluMipSize" ,"MIP size;size" ,100 ,0 ,100 );
495
496 TH2F* pM=new TH2F("RichCluMap" ,"Cluster map;x [cm];y [cm]" ,1000 ,0 ,AliHMPIDDigit::SizePcX(),1000,0,AliHMPIDDigit::SizePcY()); // maps
497 TH2F* pMipM=new TH2F("RichCluMipMap" ,"MIP map;x [cm];y [cm]" ,1000 ,0 ,AliHMPIDDigit::SizePcX(),1000,0,AliHMPIDDigit::SizePcY());
498 TH2F* pCerM=new TH2F("RichCluCerMap" ,"Ckov map;x [cm];y [cm]" ,1000 ,0 ,AliHMPIDDigit::SizePcX(),1000,0,AliHMPIDDigit::SizePcY());
499
500
501
502 for(Int_t iEvt=0;iEvt<iNevt;iEvt++){
503 pAL->GetEvent(iEvt);
504 pRL->TreeD()->GetEntry(0);
505 TClonesArray *pCluLst=new TClonesArray("AliHMPIDCluster");//tmp list of clusters for this event
506
507 for(Int_t iCh=0;iCh<7;iCh++) Dig2Clu(pRich->DigLst(iCh),pCluLst,kFALSE);//cluster finder for all chamber if any digits present
508
509 for(Int_t iClu=0;iClu<pCluLst->GetEntriesFast();iClu++){
510 AliHMPIDCluster *pClu = (AliHMPIDCluster*)pCluLst->At(iClu);
511 Int_t cfm=0; for(Int_t iDig=0;iDig<pClu->Size();iDig++) cfm+=pClu->Dig(iDig)->Ch(); //collect ckov-fee-mip structure of current cluster ?????
512 Int_t iNckov=cfm/1000000; Int_t iNfee =cfm%1000000/1000; Int_t iNmip =cfm%1000000%1000;
513
514 pQ ->Fill(pClu->Q()) ; pS ->Fill(pClu->Size()) ; pM ->Fill(pClu->X(),pClu->Y()); //all clusters
515 if(iNckov!=0 && iNfee==0 && iNmip==0) {pCerQ->Fill(pClu->Q()) ; pCerS->Fill(pClu->Size()) ; pCerM ->Fill(pClu->X(),pClu->Y());}//ckov only cluster
516 if(iNckov==0 && iNfee==0 && iNmip!=0) {pMipQ->Fill(pClu->Q()) ; pMipS->Fill(pClu->Size()) ; pMipM ->Fill(pClu->X(),pClu->Y());}//mip only cluster
517
518 }//clusters loop
519 pCluLst->Clear();delete pCluLst;
520 }//events loop
521
522 pRL->UnloadDigits(); pAL->UnloadKinematics(); pAL->UnloadHeader();
523 TCanvas *pC=new TCanvas("RichCluQA",Form("QA for cluster from %i events",iNevt),1000,900); pC->Divide(3,3);
524 pC->cd(1); pM->Draw(); pC->cd(2); pQ->Draw(); pC->cd(3); pS->Draw();
525 pC->cd(4); pMipM->Draw(); pC->cd(5); pMipQ->Draw(); pC->cd(6); pMipS->Draw();
526 pC->cd(7); pCerM->Draw(); pC->cd(8); pCerQ->Draw(); pC->cd(9); pCerS->Draw();
527}//CluQA()
528
529
530
531void AliHMPID::OccupancyPrint(Int_t iEvtNreq)
532{
533//prints occupancy for each chamber in a given event
534 Int_t iEvtNmin,iEvtNmax;
535 if(iEvtNreq==-1){
536 iEvtNmin=0;
537 iEvtNmax=gAlice->GetEventsPerRun();
538 } else {
539 iEvtNmin=iEvtNreq;iEvtNmax=iEvtNreq+1;
540 }
541
542 if(GetLoader()->GetRunLoader()->LoadHeader()) return;
543 if(GetLoader()->GetRunLoader()->LoadKinematics()) return;
544
545// Info("Occupancy","for event %i",iEvtN);
546 if(GetLoader()->LoadHits()) return;
547 if(GetLoader()->LoadDigits()) return;
548
549
550 for(Int_t iEvtN=iEvtNmin;iEvtN<iEvtNmax;iEvtN++){
551 Int_t nDigCh[7]={0,0,0,0,0,0,0};
552 Int_t iChHits[7]={0,0,0,0,0,0,0};
553 Int_t nPrim[7]={0,0,0,0,0,0,0};
554 Int_t nSec[7]={0,0,0,0,0,0,0};
555 AliInfo(Form("events processed %i",iEvtN));
556 if(GetLoader()->GetRunLoader()->GetEvent(iEvtN)) return;
557 AliStack *pStack = GetLoader()->GetRunLoader()->Stack();
558 for(Int_t iPrimN=0;iPrimN<GetLoader()->TreeH()->GetEntries();iPrimN++){//prims loop
559 GetLoader()->TreeH()->GetEntry(iPrimN);
560 for(Int_t iHitN=0;iHitN<Hits()->GetEntries();iHitN++){
561 AliHMPIDHit *pHit = (AliHMPIDHit*)Hits()->At(iHitN);
562 if(pHit->E()>0){
563 iChHits[pHit->Ch()]++;
564 if(pStack->Particle(pHit->GetTrack())->Rho()<0.01) nPrim[pHit->Ch()]++;else nSec[pHit->Ch()]++;
565 }
566 }
567 }
568
569 GetLoader()->TreeD()->GetEntry(0);
570 for(Int_t iCh=0;iCh<7;iCh++){
571 for(Int_t iDig=0;iDig<DigLst(iCh)->GetEntries();iDig++){
572 AliHMPIDDigit *pDig=(AliHMPIDDigit*)DigLst(iCh)->At(iDig);
573 nDigCh[pDig->Ch()]++;
574 }
575 Printf("Occupancy for chamber %i = %4.2f %% and charged prim tracks %i and sec. tracks %i with total %i",
576 iCh,Float_t(nDigCh[iCh])*100/AliHMPIDDigit::kPadAll,nPrim[iCh],nSec[iCh],iChHits[iCh]);
577 }
578
579
580 }//events loop
581 GetLoader()->UnloadHits();
582 GetLoader()->UnloadDigits();
583 GetLoader()->GetRunLoader()->UnloadHeader();
584 GetLoader()->GetRunLoader()->UnloadKinematics();
585}
586
587
588//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
589void AliHMPID::SummaryOfEvent(Int_t iEvtN) const
590{
591//prints a summary for a given event
592 AliInfo(Form("Summary of event %i",iEvtN));
593 GetLoader()->GetRunLoader()->GetEvent(iEvtN);
594 if(GetLoader()->GetRunLoader()->LoadHeader()) return;
595 if(GetLoader()->GetRunLoader()->LoadKinematics()) return;
596 AliStack *pStack=GetLoader()->GetRunLoader()->Stack();
597
598 AliGenEventHeader* pGenHeader = gAlice->GetHeader()->GenEventHeader();
599 if(pGenHeader->InheritsFrom("AliGenHijingEventHeader")) {
600 AliInfo(Form(" Hijing event with impact parameter b = %.2f (fm)",((AliGenHijingEventHeader*) pGenHeader)->ImpactParameter()));
601 }
602 Int_t nChargedPrimaries=0;
603 for(Int_t i=0;i<pStack->GetNtrack();i++) {
604 TParticle *pParticle = pStack->Particle(i);
605 if(pParticle->IsPrimary()&&pParticle->GetPDG()->Charge()!=0) nChargedPrimaries++;
606 }
607 AliInfo(Form("Total number of primaries %i",pStack->GetNprimary()));
608 AliInfo(Form("Total number of charged primaries %i",nChargedPrimaries));
609 AliInfo(Form("Total n. of tracks in stack(+sec) %i",pStack->GetNtrack()));
610 GetLoader()->GetRunLoader()->UnloadHeader();
611 GetLoader()->GetRunLoader()->UnloadKinematics();
612}
613//++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++