From Magali Estienne.
[u/mrichter/AliRoot.git] / EVE / EveDet / AliEveEMCALData.cxx
1 //*********************************************************************
2 // - AliEVE implementation -
3 // Fill containers for visualisation of EMCAL data structures
4 //    - read and store MC Hits
5 //    - read and store digits from esds or runloader
6 //    - read and store clusters from esds or runloader 
7 //
8 // Author: Magali Estienne (magali.estienne@cern.ch)
9 // June 30 2008
10 //*********************************************************************
11
12 #include <TObject.h>
13 #include <TEveUtil.h>
14 #include <TTree.h>
15 #include <TBranch.h>
16 #include <TObjArray.h>
17 #include <TRefArray.h>
18 #include <TClonesArray.h>
19 #include <TEvePointSet.h>
20
21 #include "AliRun.h"
22 #include "AliRunLoader.h"
23 #include "AliEMCALLoader.h"
24 #include "AliESDEvent.h"
25 #include "AliESDVertex.h"
26 #include "AliEMCAL.h"
27 #include "AliEMCALGeometry.h"
28 #include "AliEMCALHit.h"
29 #include "AliEMCALDigit.h"
30 #include "AliEMCALRecPoint.h"
31 #include "AliESDCaloCells.h"
32 #include "AliESDCaloCluster.h"
33
34 #include "AliEveEMCALData.h"
35 #include "AliEveEMCALSModule.h"
36 #include "AliEveEMCALSModuleData.h"
37
38 ClassImp(AliEveEMCALData)
39
40 //______________________________________________________________________________
41 AliEveEMCALData::AliEveEMCALData():
42   TObject(),
43   TEveRefCnt(),
44   fEmcal(0x0),
45   fGeom(0x0),
46   fNode(0x0),
47   fHMatrix(0),
48   fTree(0x0),
49   fESD(0x0),
50   fNsm(12),
51   fNsmfull(10),
52   fNsmhalf(2),
53   fSM(12),
54   fSMfull(10),
55   fSMhalf(2),
56   fRunLoader(0),
57   fDebug(0),
58   fPoint(0)
59 {
60   
61   //
62   // Constructor
63   //
64   //   for(Int_t i=0; i<12; i++)
65   //     fSM[i] = 0x0;
66   //   for(Int_t i=0; i<10; i++)
67   //     fSMfull[i] = 0x0;
68   //   fSMhalf[0] = fSMhalf[1] = 0x0;
69   //  InitEMCALGeom();
70   CreateAllSModules();
71
72
73 }
74
75 //______________________________________________________________________________
76 AliEveEMCALData::AliEveEMCALData(AliRunLoader* rl, TGeoNode* node, TGeoHMatrix* m):
77   TObject(),
78   TEveRefCnt(),
79   fEmcal(0x0),
80   fGeom(0x0),
81   fNode(node),
82   fHMatrix(m),
83   fTree(0x0),
84   fESD(0x0),
85   fNsm(12),
86   fNsmfull(10),
87   fNsmhalf(2),
88   fSM(12),
89   fSMfull(10),
90   fSMhalf(2),
91   fRunLoader(rl),
92   fDebug(0),
93   fPoint(0)
94 {
95   
96   //
97   // Constructor
98   //
99 //   for(Int_t i=0; i<12; i++)
100 //     fSM[i] = 0x0;
101 //   for(Int_t i=0; i<10; i++)
102 //     fSMfull[i] = 0x0;
103 //   fSMhalf[0] = fSMhalf[1] = 0x0;
104   InitEMCALGeom(rl);
105   CreateAllSModules();
106
107
108 }
109
110 //______________________________________________________________________________
111 AliEveEMCALData::~AliEveEMCALData()
112 {
113   //
114   // Destructor
115   //
116
117   DeleteSuperModules();
118   delete fTree;
119   delete fEmcal;
120   delete fGeom;
121   delete fNode;
122   delete fHMatrix;
123   delete fPoint;
124 }
125
126 //______________________________________________________________________________
127 void AliEveEMCALData::Reset()
128 {
129
130 }
131
132 //______________________________________________________________________________
133 AliEveEMCALData::AliEveEMCALData(const AliEveEMCALData &edata) :
134   TObject(edata),
135   TEveRefCnt(edata),
136   fEmcal(edata.fEmcal),
137   fGeom(edata.fGeom),
138   fNode(edata.fNode),
139   fHMatrix(edata.fHMatrix),
140   fTree(edata.fTree),
141   fESD(edata.fESD),
142   fNsm(edata.fNsm),
143   fNsmfull(edata.fNsmfull),
144   fNsmhalf(edata.fNsmhalf),
145   fSM(edata.fSM),
146   fSMfull(edata.fSMfull),
147   fSMhalf(edata.fSMhalf),
148   fRunLoader(edata.fRunLoader),
149   fDebug(edata.fDebug),
150   fPoint(edata.fPoint)
151 {
152   //
153   // Copy constructor
154   //
155   InitEMCALGeom(edata.fRunLoader);
156   CreateAllSModules();
157 }
158
159 //______________________________________________________________________________
160 AliEveEMCALData& AliEveEMCALData::operator=(const AliEveEMCALData &edata)
161 {
162   //
163   // Assignment operator
164   //
165
166   if (this != &edata) {
167
168   }
169
170   return *this;
171
172 }
173
174 //______________________________________________________________________________
175 void AliEveEMCALData::SetTree(TTree* tree)
176 {
177
178   // Set digit-tree to be used for digit retrieval. Data is loaded on
179   // demand.
180
181   fTree = tree;
182
183 }
184
185 //______________________________________________________________________________
186 void AliEveEMCALData::SetESD(AliESDEvent* esd)
187 {
188   fESD = esd;
189 }
190
191 //______________________________________________________________________________
192 void AliEveEMCALData::SetNode(TGeoNode* node)
193 {
194   fNode = node;
195 }
196
197 //______________________________________________________________________________
198 void AliEveEMCALData::InitEMCALGeom(AliRunLoader* rl)
199 {
200 //   // Handle an individual point selection from GL.
201
202 //   fParent->SpawnEditor();
203   fEmcal = (AliEMCAL*) rl->GetAliRun()->GetDetector("EMCAL");
204   fGeom  = (AliEMCALGeometry*) fEmcal->GetGeometry();
205
206 }
207
208 //______________________________________________________________________________
209 void AliEveEMCALData::GetGeomInfo(Int_t id, Int_t &iSupMod, Double_t& x, Double_t& y, Double_t& z)
210 {
211   // Get geometrical information from hit/digit/cluster absolute id
212
213   Int_t iTower  =  0 ;
214   Int_t iIphi   =  0 ;
215   Int_t iIeta   =  0 ;
216   Int_t iphi    =  0 ;
217   Int_t ieta    =  0 ;
218
219   //Geometry methods
220   fGeom->GetCellIndex(id,iSupMod,iTower,iIphi,iIeta);
221   //Gives SuperModule and Tower numbers
222   fGeom->RelPosCellInSModule(id, x, y, z);
223
224 }
225
226 //______________________________________________________________________________
227 void  AliEveEMCALData::CreateAllSModules()
228 {
229
230   // Create all fNsm super modules
231   for(Int_t sm = 0; sm < fNsm; sm++)
232     CreateSModule(sm);
233
234 }
235
236 //______________________________________________________________________________
237 void  AliEveEMCALData::CreateSModule(Int_t sm)
238 {
239   // Create super-module-data for sm if it does not exist already.
240   if(fSM[sm] == 0) fSM[sm] = new AliEveEMCALSModuleData(sm,fGeom,fNode,fHMatrix);
241   if(fSMfull[sm] == 0 && sm < 10) fSMfull[sm] = new AliEveEMCALSModuleData(sm,fGeom,fNode,fHMatrix);
242   if(fSMhalf[sm-10] == 0 && sm > 10) fSMhalf[sm-10] = new AliEveEMCALSModuleData(sm,fGeom,fNode,fHMatrix);
243 }
244
245 //______________________________________________________________________________
246 void AliEveEMCALData::DropAllSModules()
247 {
248   // Drop data of all existing sectors.
249
250   for (Int_t sm = 0; sm < fNsm; sm++) {
251     if (fSM[sm] != 0)
252       fSM[sm]->DropData();
253   }
254 }
255
256 //______________________________________________________________________________
257 void AliEveEMCALData::DeleteSuperModules()
258 {
259   //
260   // delete all super module data
261   //
262
263   for (Int_t sm = 0; sm < fNsm; sm++)
264     {
265       fSM[sm] = 0;
266       delete fSM[sm];
267     }
268
269   for(Int_t smf = 0; smf < fNsmfull; smf++) 
270     {
271       fSMfull[smf] = 0;
272       delete fSMfull[smf];
273     }
274
275   for(Int_t smh = 0; smh < fNsmhalf; smh++)
276     {
277       fSMhalf[smh] = 0;
278       delete fSMhalf[smh];
279     }
280   
281 }
282
283 //______________________________________________________________________________
284 void AliEveEMCALData::LoadHits(TTree* t)
285 {
286
287   // With TTreeH
288
289   // These are global coordinates !
290   char form[1000];
291   const char *selection = "";
292   const char *varexp = "fX:fY:fZ";
293   sprintf(form,"EMCAL Hits '%s'", selection);
294   fPoint = new TEvePointSet(form);
295
296   TEvePointSelector ps(t, fPoint, varexp, selection);
297   ps.Select();
298
299   if (fPoint->Size() == 0) {
300     Warning("emcal_hits", Form("No hits match '%s'", selection));
301     delete fPoint;
302     //    return 0;
303   }
304
305
306 //   TObjArray *harr=NULL;
307 //   TBranch *hbranch=t->GetBranch("EMCAL");
308 //   hbranch->SetAddress(&harr);
309   
310 //   if(hbranch->GetEvent(0)) {
311 //     for(Int_t ih = 0; ih < harr->GetEntriesFast(); ih++) {
312 //       AliEMCALHit* hit =(AliEMCALHit*)harr->UncheckedAt(ih);
313 //       if(hit != 0){
314 //      cout << "Hit info " << hit->GetId() << " " << hit->GetEnergy() << endl;
315 //      fSM[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
316 //       }
317 //     }
318 //   }
319 //   //********************************
320 //   // To be completed !!!
321 //   //********************************
322
323 }
324
325 //______________________________________________________________________________
326 void AliEveEMCALData::LoadHitsFromEMCALLoader(AliEMCALLoader* emcl)
327 {
328
329   // Need to know the local position in each sm
330
331   // With EMCALLoader
332   AliEMCALHit* hit;
333   
334   //Fill array of hits                                                                        
335   TClonesArray *hits = (TClonesArray*)emcl->Hits();
336
337   //Get hits from the list                                                                    
338   for(Int_t ihit = 0; ihit< hits->GetEntries();ihit++){
339     //cout<<">> idig "<<idig<<endl;                                                             
340     hit = static_cast<AliEMCALHit *>(hits->At(ihit)) ;
341     
342     if(hit != 0){
343       cout << "Hit info " << hit->GetId() << " " << hit->GetEnergy() << endl;
344     }
345   }
346   
347   //********************************
348   // To be completed !!!
349   //********************************
350
351
352
353 }
354
355 //______________________________________________________________________________
356 void AliEveEMCALData::LoadDigits(TTree *t)
357 {
358   //
359   // load digits from the TreeD
360   //
361
362   TClonesArray *digits = 0;
363   t->SetBranchAddress("EMCAL", &digits);
364   t->GetEntry(0);
365
366   Int_t nEnt = digits->GetEntriesFast();
367   cout << "nEnt: " << nEnt << endl;
368   AliEMCALDigit * dig;
369
370   Float_t amp   = -1 ;
371   Int_t id      = -1 ;
372   Int_t iSupMod =  0 ;
373   Double_t x, y, z;
374
375   for (Int_t idig = 0; idig < nEnt; idig++)
376     {
377       dig = static_cast<AliEMCALDigit *>(digits->At(idig));
378       
379       if(dig != 0) {
380         id   = dig->GetId() ; //cell (digit) label
381         amp  = dig->GetAmp(); //amplitude in cell (digit)
382         
383         GetGeomInfo(id,iSupMod,x,y,z);
384
385         fSM[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
386 //      fSM[iSupMod]->SaveDigit(dig);
387 //      if(iSupMod<fNsmfull) fSMfull[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
388 //      if(iSupMod>fNsmfull) fSMhalf[iSupMod-10]->RegisterDigit(id,iSupMod,amp,x,y,z);
389       }
390       else {
391         cout << "Digit object empty" << endl;
392         return;
393       }
394     } // end loop digits
395   cout << "after loop on digits !" << endl;
396 }
397
398 //______________________________________________________________________________
399 void AliEveEMCALData::LoadDigitsFromEMCALLoader(AliEMCALLoader* emcl)
400 {
401
402   //
403   // load digits from EMCAL Loader
404   //
405
406   AliEMCALDigit* dig;
407   
408   //Fill array of digits                                                                        
409   TClonesArray *digits = (TClonesArray*)emcl->Digits();
410   
411   //Get digits from the list  
412   
413   Float_t amp   = -1 ;
414   Int_t id      = -1 ;
415   Int_t iSupMod =  0 ;
416   Double_t x, y, z;
417                                                
418    for(Int_t idig = 0; idig< digits->GetEntries();idig++){
419
420      dig = static_cast<AliEMCALDigit *>(digits->At(idig)) ;
421
422      if(dig != 0){
423        if(fDebug>1) cout << "Digit info " << dig->GetId() << " " << dig->GetAmp() << endl;
424        id   = dig->GetId() ; //cell (digit) label
425        amp  = dig->GetAmp(); //amplitude in cell (digit)
426        
427        GetGeomInfo(id,iSupMod,x,y,z);
428        
429        fSM[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
430      }
431       else {
432         cout << "Digit object empty" << endl;
433         return;
434       }
435    } // end loop on digits
436    
437 }
438
439 //______________________________________________________________________________
440 void AliEveEMCALData::LoadDigitsFromESD()
441 {
442   //
443   // Get digit information from esd
444   //
445   
446   AliESDCaloCells &cells= *(fESD->GetEMCALCells());
447   Int_t ncell = cells.GetNumberOfCells() ;  
448   Int_t type = cells.GetType();
449   Float_t pos[3] ; 
450   Int_t iSupMod =  0 ;
451   Double_t x, y, z;
452
453   // Extract digit information from the ESDs
454   for (Int_t icell=  0; icell <  ncell; icell++) 
455     {
456       Int_t id   = cells.GetCellNumber(icell);
457       Float_t amp  = cells.GetAmplitude(icell);
458
459       GetGeomInfo(id,iSupMod,x,y,z);
460
461       fSM[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
462       if(iSupMod<fNsmfull) fSMfull[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
463       if(iSupMod>fNsmfull) fSMhalf[iSupMod-10]->RegisterDigit(id,iSupMod,amp,x,y,z);
464
465     } // end loop cells
466 }
467
468 //______________________________________________________________________________
469 void AliEveEMCALData::LoadRecPoints(TTree* t)
470 {
471   //*************************************************
472   // To be improved !!!!!
473   // Size and shape of cluster to be implemented
474   // 
475   //*************************************************
476         
477   // From TTreeR
478   TObjArray *carr=NULL;
479   TBranch *cbranch=t->GetBranch("EMCALECARP");
480   cbranch->SetAddress(&carr);
481   
482   if(cbranch->GetEvent(0)) {
483     for(Int_t ic = 0; ic < carr->GetEntriesFast(); ic++) {
484       AliEMCALRecPoint* rp =(AliEMCALRecPoint*)carr->UncheckedAt(ic);
485       if(rp){
486         cout << "RecPoint info " << rp->GetAbsId() << " " << rp->GetEnergy() << endl;
487         Int_t iSupMod = rp->GetSuperModuleNumber();
488         Float_t amp = rp->GetEnergy();
489         TVector3 lpos;
490         rp->GetLocalPosition(lpos);
491
492         fSM[iSupMod]->RegisterCluster(iSupMod,amp,lpos[0],lpos[1],lpos[2]);
493       }
494     }
495   }
496   
497 }
498
499 //______________________________________________________________________________
500 void AliEveEMCALData::LoadRecPointsFromEMCALLoader(AliEMCALLoader* emcl)
501 {
502   //*************************************************
503   // To be improved !!!!!
504   // Size and shape of cluster to be implemented
505   // 
506   //*************************************************
507
508   // From EMCALLoader
509   AliEMCALRecPoint* rp;
510   
511   //Fill array of clusters                                                                        
512   TClonesArray *clusters = (TClonesArray*)emcl->RecPoints();
513   
514   //Get clusters from the list                                                                    
515   for(Int_t iclu = 0; iclu< clusters->GetEntries();iclu++){
516
517     rp = static_cast<AliEMCALRecPoint *>(clusters->At(iclu)) ;
518     
519     if(rp){
520        cout << "RecPoint info " << rp->GetAbsId() << " " << rp->GetEnergy() << endl;
521        Int_t iSupMod = rp->GetSuperModuleNumber();
522        Float_t amp = rp->GetEnergy();
523        TVector3 lpos;
524        rp->GetLocalPosition(lpos);
525        
526        fSM[iSupMod]->RegisterCluster(iSupMod,amp,lpos[0],lpos[1],lpos[2]);
527     }
528   }
529   
530 }
531
532 //______________________________________________________________________________
533 void AliEveEMCALData::LoadRecPointsFromESD()
534 {
535   Int_t iSupMod =  0 ;
536   Double_t x, y, z;
537   Int_t iSM =  0 ;
538   Int_t iT  =  0 ;
539   Int_t iIp   =  0 ;
540   Int_t iIe   =  0 ;
541   Int_t ip    =  0 ;
542   Int_t ie    =  0 ;
543   Double_t xd, yd, zd;
544   Float_t pos[3] ; 
545   
546   // Get reconstructed vertex position
547   AliESDVertex* primVertex =(AliESDVertex*) fESD->GetVertex();
548   Double_t vertex_position[3] ; 
549   primVertex->GetXYZ(vertex_position) ; 
550
551   //Get the CaloClusters
552   //select EMCAL clusters only 
553   TRefArray * caloClusters  = new TRefArray();
554   fESD->GetEMCALClusters(caloClusters);
555   Int_t nclus = caloClusters->GetEntries();
556   cout << "nclus: " << nclus << endl; 
557   
558   if(!caloClusters) return;
559
560   for (Int_t iclus =  0; iclus <  nclus; iclus++) 
561     {
562       AliESDCaloCluster *clus = (AliESDCaloCluster *) caloClusters->At(iclus) ; 
563       //Get the cluster info
564       
565       Float_t energy = clus->E() ;
566       Int_t   eneInt = (Int_t) energy*500+0.5;
567       Float_t disp   = clus->GetClusterDisp() ;
568       Int_t   iprim  = clus->GetLabel();
569       
570       clus->GetPosition(pos) ; // Global position
571       TVector3 vpos(pos[0],pos[1],pos[2]) ;
572       TLorentzVector p4 ;
573       TVector3 p3;
574       clus->GetMomentum(p4,vertex_position);
575       p3.SetXYZ(p4[0],p4[1],p4[2]);
576       Float_t eta = p3.Eta();
577       Float_t phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi());
578       
579       Int_t mult = clus->GetNCells() ;
580       if(fDebug>1) {
581         cout << "In cluster: " << iclus << ", ncells: " << mult << ", energy: " << 
582           eneInt << ", disp: " << disp << endl;
583         cout << "Cluster " << iclus << ", eta: " << eta << ", phi: " << phi << endl;
584       }
585
586       Int_t clusId = 0;
587       fGeom->GetAbsCellIdFromEtaPhi(eta,phi,clusId);
588       if(fDebug>1) {
589         cout << "Abs Cluster Id: " << clusId << ", xc: " << pos[0] << 
590           ", yc: " << pos[1] << ", zc: " << pos[2] << endl;
591       }
592
593       GetGeomInfo(clusId,iSupMod,x,y,z);
594       
595       //******** Not used yet but will come  ********
596       Int_t     digMult = clus->GetNCells() ;
597       UShort_t *digID   = clus->GetCellsAbsId() ;
598       for(Int_t i=0; i<digMult; i++){
599         fGeom->RelPosCellInSModule(digID[i], xd, yd, zd);
600         //Geometry methods
601         fGeom->GetCellIndex(digID[i],iSM,iT,iIp,iIe);
602         //Gives SuperModule and Tower numbers
603
604       } // end digit loop
605       //*********************************************
606
607       fSM[iSupMod]->RegisterCluster(iSM,energy,x,y,z);
608
609     } // end cluster loop
610 }
611
612 //______________________________________________________________________________
613 AliEveEMCALSModuleData* AliEveEMCALData::GetSModuleData(Int_t sm)
614 {
615   //
616   // return super module data
617   //
618
619   if (sm < 0 || sm > fNsm) 
620     {
621       printf("The number of super modules must be lower or equal to %d",fNsm);
622       return 0;
623     }
624   
625   return fSM[sm];
626 }
627
628 //______________________________________________________________________________
629 void AliEveEMCALData::LoadRaw()
630 {
631
632 }
633