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
8 // Author: Magali Estienne (magali.estienne@cern.ch)
10 //*********************************************************************
16 #include <TObjArray.h>
17 #include <TRefArray.h>
18 #include <TClonesArray.h>
19 #include <TEvePointSet.h>
22 #include "AliRunLoader.h"
23 #include "AliEMCALLoader.h"
24 #include "AliESDEvent.h"
25 #include "AliESDVertex.h"
27 #include "AliEMCALGeometry.h"
28 #include "AliEMCALHit.h"
29 #include "AliEMCALDigit.h"
31 #include "AliEMCALRecPoint.h"
32 #include "AliESDCaloCells.h"
33 #include "AliESDCaloCluster.h"
35 #include "AliEveEMCALData.h"
36 #include "AliEveEMCALSModule.h"
37 #include "AliEveEMCALSModuleData.h"
39 ClassImp(AliEveEMCALData)
41 //______________________________________________________________________________
42 AliEveEMCALData::AliEveEMCALData():
65 // for(Int_t i=0; i<12; i++)
67 // for(Int_t i=0; i<10; i++)
69 // fSMhalf[0] = fSMhalf[1] = 0x0;
76 //______________________________________________________________________________
77 AliEveEMCALData::AliEveEMCALData(AliRunLoader* rl, TGeoNode* node, TGeoHMatrix* m):
100 // for(Int_t i=0; i<12; i++)
102 // for(Int_t i=0; i<10; i++)
104 // fSMhalf[0] = fSMhalf[1] = 0x0;
111 //______________________________________________________________________________
112 AliEveEMCALData::~AliEveEMCALData()
118 DeleteSuperModules();
120 // delete fEmcal; // deleted by run-loader
127 //______________________________________________________________________________
128 void AliEveEMCALData::Reset()
133 //______________________________________________________________________________
134 AliEveEMCALData::AliEveEMCALData(const AliEveEMCALData &edata) :
137 fEmcal(edata.fEmcal),
140 fHMatrix(edata.fHMatrix),
144 fNsmfull(edata.fNsmfull),
145 fNsmhalf(edata.fNsmhalf),
147 fSMfull(edata.fSMfull),
148 fSMhalf(edata.fSMhalf),
149 fRunLoader(edata.fRunLoader),
150 fDebug(edata.fDebug),
156 InitEMCALGeom(edata.fRunLoader);
160 //______________________________________________________________________________
161 AliEveEMCALData& AliEveEMCALData::operator=(const AliEveEMCALData &edata)
164 // Assignment operator
167 if (this != &edata) {
175 //______________________________________________________________________________
176 void AliEveEMCALData::SetTree(TTree* tree)
179 // Set digit-tree to be used for digit retrieval.
180 // Data is loaded on demand.
187 //______________________________________________________________________________
188 void AliEveEMCALData::SetESD(AliESDEvent* esd)
197 //______________________________________________________________________________
198 void AliEveEMCALData::SetNode(TGeoNode* node)
207 //______________________________________________________________________________
208 void AliEveEMCALData::InitEMCALGeom(AliRunLoader* rl)
211 // Set data members for EMCAL geometry
214 fEmcal = (AliEMCAL*) rl->GetAliRun()->GetDetector("EMCAL");
215 fGeom = (AliEMCALGeometry*) fEmcal->GetGeometry();
219 //______________________________________________________________________________
220 void AliEveEMCALData::GetGeomInfo(Int_t id, Int_t &iSupMod, Double_t& x, Double_t& y, Double_t& z)
223 // Get geometrical information from hit/digit/cluster absolute id
231 fGeom->GetCellIndex(id,iSupMod,iTower,iIphi,iIeta);
232 //Gives SuperModule and Tower numbers
233 fGeom->RelPosCellInSModule(id, x, y, z);
237 //______________________________________________________________________________
238 void AliEveEMCALData::CreateAllSModules()
241 // Create all fNsm super modules
244 for(Int_t sm = 0; sm < fNsm; sm++)
249 //______________________________________________________________________________
250 void AliEveEMCALData::CreateSModule(Int_t sm)
253 // Create super-module-data for SM if it does not exist already.
256 if(fSM[sm] == 0) fSM[sm] = new AliEveEMCALSModuleData(sm,fGeom,fNode,fHMatrix);
257 if(fSMfull[sm] == 0 && sm < 10) fSMfull[sm] = new AliEveEMCALSModuleData(sm,fGeom,fNode,fHMatrix);
258 if(fSMhalf[sm-10] == 0 && sm > 10) fSMhalf[sm-10] = new AliEveEMCALSModuleData(sm,fGeom,fNode,fHMatrix);
261 //______________________________________________________________________________
262 void AliEveEMCALData::DropAllSModules()
265 // Drop data of all existing sectors.
268 for (Int_t sm = 0; sm < fNsm; sm++) {
274 //______________________________________________________________________________
275 void AliEveEMCALData::DeleteSuperModules()
278 // Delete all super module data
281 for (Int_t sm = 0; sm < fNsm; sm++)
287 for(Int_t smf = 0; smf < fNsmfull; smf++)
293 for(Int_t smh = 0; smh < fNsmhalf; smh++)
301 //______________________________________________________________________________
302 void AliEveEMCALData::LoadHits(TTree* t)
305 // Get hit information from RunLoader
309 // These are global coordinates !
311 const char *selection = "";
312 const char *varexp = "fX:fY:fZ";
313 sprintf(form,"EMCAL Hits '%s'", selection);
314 fPoint = new TEvePointSet(form);
316 TEvePointSelector ps(t, fPoint, varexp, selection);
319 if (fPoint->Size() == 0) {
320 Warning("emcal_hits", Form("No hits match '%s'", selection));
326 TObjArray *harr=NULL;
327 TBranch *hbranch=t->GetBranch("EMCAL");
328 hbranch->SetAddress(&harr);
330 if(hbranch->GetEvent(0)) {
331 for(Int_t ih = 0; ih < harr->GetEntriesFast(); ih++) {
332 AliEMCALHit* hit =(AliEMCALHit*)harr->UncheckedAt(ih);
334 if(fDebug>1) cout << "Hit info " << hit->GetId() << " " << hit->GetEnergy() << endl;
335 Int_t id = hit->GetId();
336 // These are local coordinates
337 Double_t xl = 0.; Double_t yl = 0.; Double_t zl = 0.;
338 // Get global coordinates
339 Double_t x = hit->X();
340 Double_t y = hit->Y();
341 Double_t z = hit->Z();
342 Double_t amp = hit->GetEnergy();
345 GetGeomInfo(id,iSupMod,xl,yl,zl);
346 fSM[iSupMod]->RegisterHit(id,iSupMod,amp,x,y,z);
352 //______________________________________________________________________________
353 void AliEveEMCALData::LoadHitsFromEMCALLoader(AliEMCALLoader* emcl)
356 // Get hit information from EMCAL Loader
362 TClonesArray *hits = (TClonesArray*)emcl->Hits();
364 //Get hits from the list
365 for(Int_t ihit = 0; ihit< hits->GetEntries();ihit++){
367 hit = static_cast<AliEMCALHit *>(hits->At(ihit)) ;
370 if(fDebug>1) cout << "Hit info " << hit->GetId() << " " << hit->GetEnergy() << endl;
372 Int_t id = hit->GetId();
373 // These are local coordinates
374 Double_t xl = 0.; Double_t yl = 0.; Double_t zl = 0.;
375 // Get global coordinates
376 Double_t x = hit->X();
377 Double_t y = hit->Y();
378 Double_t z = hit->Z();
379 Double_t amp = hit->GetEnergy();
382 GetGeomInfo(id,iSupMod,xl,yl,zl);
383 fSM[iSupMod]->RegisterHit(id,iSupMod,amp,x,y,z);
389 //______________________________________________________________________________
390 void AliEveEMCALData::LoadDigits(TTree *t)
393 // Get digit information from RunLoader
396 TClonesArray *digits = 0;
397 t->SetBranchAddress("EMCAL", &digits);
400 Int_t nEnt = digits->GetEntriesFast();
401 cout << "nEnt: " << nEnt << endl;
404 // Double_t amp = -1 ;
405 Double_t ampInt = -1 ;
410 for (Int_t idig = 0; idig < nEnt; idig++)
412 dig = static_cast<AliEMCALDigit *>(digits->At(idig));
415 id = dig->GetId() ; //cell (digit) label
417 ampInt = dig->GetAmp(); //amplitude in cell (digit)
419 // amp = ampInt*0.0153; // To be modified with correct OCDB conversion
421 GetGeomInfo(id,iSupMod,x,y,z);
424 // fSM[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
425 // // fSM[iSupMod]->SaveDigit(dig);
426 // // if(iSupMod<fNsmfull) fSMfull[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
427 // // if(iSupMod>fNsmfull) fSMhalf[iSupMod-10]->RegisterDigit(id,iSupMod,amp,x,y,z);
428 fSM[iSupMod]->RegisterDigit(id,iSupMod,ampInt,x,y,z);
429 // fSM[iSupMod]->SaveDigit(dig);
430 // if(iSupMod<fNsmfull) fSMfull[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
431 // if(iSupMod>fNsmfull) fSMhalf[iSupMod-10]->RegisterDigit(id,iSupMod,amp,x,y,z);
434 cout << "Digit object empty" << endl;
438 cout << "after loop on digits !" << endl;
441 //______________________________________________________________________________
442 void AliEveEMCALData::LoadDigitsFromEMCALLoader(AliEMCALLoader* emcl)
446 // Get digit information from EMCAL Loader
451 //Fill array of digits
452 TClonesArray *digits = (TClonesArray*)emcl->Digits();
454 //Get digits from the list
456 // Double_t amp = -1 ;
457 Double_t ampInt = -1 ;
462 for(Int_t idig = 0; idig< digits->GetEntries();idig++){
464 dig = static_cast<AliEMCALDigit *>(digits->At(idig)) ;
467 if(fDebug>1) cout << "Digit info " << dig->GetId() << " " << dig->GetAmp() << endl;
468 id = dig->GetId() ; //cell (digit) label
470 ampInt = dig->GetAmp(); //amplitude in cell (digit)
472 // amp = ampInt*0.0153.; // To be modified with correct OCDB conversion
474 GetGeomInfo(id,iSupMod,x,y,z);
477 // fSM[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
479 fSM[iSupMod]->RegisterDigit(id,iSupMod,ampInt,x,y,z);
482 cout << "Digit object empty" << endl;
485 } // end loop on digits
489 //______________________________________________________________________________
490 void AliEveEMCALData::LoadDigitsFromESD()
493 // Get digit information from esd
496 AliESDCaloCells &cells= *(fESD->GetEMCALCells());
497 Int_t ncell = cells.GetNumberOfCells() ;
501 // Extract digit information from the ESDs
502 for (Int_t icell= 0; icell < ncell; icell++)
504 Int_t id = cells.GetCellNumber(icell);
506 Double_t ampInt = cells.GetAmplitude(icell);
508 // Double_t amp = ampInt*0.0153; // To be modified with correct OCDB conversion
510 GetGeomInfo(id,iSupMod,x,y,z);
513 // fSM[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
514 // if(iSupMod<fNsmfull) fSMfull[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
515 // if(iSupMod>fNsmfull) fSMhalf[iSupMod-10]->RegisterDigit(id,iSupMod,amp,x,y,z);
517 fSM[iSupMod]->RegisterDigit(id,iSupMod,ampInt,x,y,z);
518 if(iSupMod<fNsmfull) fSMfull[iSupMod]->RegisterDigit(id,iSupMod,ampInt,x,y,z);
519 if(iSupMod>fNsmfull) fSMhalf[iSupMod-10]->RegisterDigit(id,iSupMod,ampInt,x,y,z);
524 //______________________________________________________________________________
525 void AliEveEMCALData::LoadRecPoints(TTree* t)
528 // Get rec point information from RunLoader
531 //*************************************************
532 // To be improved !!!!!
533 // Size and shape of cluster to be implemented
535 //*************************************************
538 TObjArray *carr=NULL;
539 TBranch *cbranch=t->GetBranch("EMCALECARP");
540 cbranch->SetAddress(&carr);
542 if(cbranch->GetEvent(0)) {
543 for(Int_t ic = 0; ic < carr->GetEntriesFast(); ic++) {
544 AliEMCALRecPoint* rp =(AliEMCALRecPoint*)carr->UncheckedAt(ic);
546 if(fDebug>1) cout << "RecPoint info " << rp->GetAbsId() << " " << rp->GetEnergy() << endl;
547 Int_t iSupMod = rp->GetSuperModuleNumber();
549 Double_t amp = (Double_t)rp->GetEnergy();
551 Double_t ampInt = amp/0.0153; // To be modified with correct OCDB conversion
553 rp->GetLocalPosition(lpos);
556 // fSM[iSupMod]->RegisterCluster(iSupMod,amp,lpos[0],lpos[1],lpos[2]);
558 fSM[iSupMod]->RegisterCluster(iSupMod,ampInt,lpos[0],lpos[1],lpos[2]);
565 //______________________________________________________________________________
566 void AliEveEMCALData::LoadRecPointsFromEMCALLoader(AliEMCALLoader* emcl)
569 // Get rec point information from EMCAL Loader
572 //*************************************************
573 // To be improved !!!!!
574 // Size and shape of cluster to be implemented
576 //*************************************************
579 AliEMCALRecPoint* rp;
581 //Fill array of clusters
582 TClonesArray *clusters = (TClonesArray*)emcl->RecPoints();
584 //Get clusters from the list
585 for(Int_t iclu = 0; iclu< clusters->GetEntries();iclu++){
587 rp = static_cast<AliEMCALRecPoint *>(clusters->At(iclu)) ;
590 if(fDebug>1) cout << "RecPoint info " << rp->GetAbsId() << " " << rp->GetEnergy() << endl;
591 Int_t iSupMod = rp->GetSuperModuleNumber();
592 Double_t amp = (Double_t)rp->GetEnergy();
593 Double_t ampInt = amp/0.0153; // To be modified with correct OCDB conversion
595 rp->GetLocalPosition(lpos);
598 // fSM[iSupMod]->RegisterCluster(iSupMod,amp,lpos[0],lpos[1],lpos[2]);
600 fSM[iSupMod]->RegisterCluster(iSupMod,ampInt,lpos[0],lpos[1],lpos[2]);
606 //______________________________________________________________________________
607 void AliEveEMCALData::LoadRecPointsFromESD()
610 // Get cluster information from esd
622 // Get reconstructed vertex position
623 AliESDVertex* primVertex =(AliESDVertex*) fESD->GetVertex();
624 Double_t vertex_position[3] ;
625 primVertex->GetXYZ(vertex_position) ;
627 //Get the CaloClusters
628 //select EMCAL clusters only
629 TRefArray * caloClusters = new TRefArray();
630 fESD->GetEMCALClusters(caloClusters);
631 Int_t nclus = caloClusters->GetEntries();
632 cout << "nclus: " << nclus << endl;
634 if(!caloClusters) return;
636 for (Int_t iclus = 0; iclus < nclus; iclus++)
638 AliESDCaloCluster *clus = (AliESDCaloCluster *) caloClusters->At(iclus) ;
639 //Get the cluster info
641 Double_t energy = clus->E() ;
643 // Int_t eneInt = (Int_t)energy*500+0.5;
644 Double_t eneInt = energy/0.0153; // To be modified with correct OCDB conversion
645 Double_t disp = clus->GetClusterDisp() ;
647 clus->GetPosition(pos) ; // Global position
648 TVector3 vpos(pos[0],pos[1],pos[2]) ;
651 clus->GetMomentum(p4,vertex_position);
652 p3.SetXYZ(p4[0],p4[1],p4[2]);
653 Double_t eta = p3.Eta();
654 Double_t phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi());
656 Int_t mult = clus->GetNCells() ;
658 cout << "In cluster: " << iclus << ", ncells: " << mult << ", energy : " << energy <<
659 ", disp: " << disp << endl;
660 cout << "Cluster " << iclus << ", eta: " << eta << ", phi: " << phi << endl;
664 fGeom->GetAbsCellIdFromEtaPhi(eta,phi,clusId);
666 cout << "Abs Cluster Id: " << clusId << ", xc: " << pos[0] <<
667 ", yc: " << pos[1] << ", zc: " << pos[2] << endl;
670 GetGeomInfo(clusId,iSupMod,x,y,z);
672 //******** Not used yet but will come ********
673 // AliESDCaloCells &cells= *(fESD->GetEMCALCells());
674 Int_t digMult = clus->GetNCells() ;
675 UShort_t *digID = clus->GetCellsAbsId() ;
676 for(Int_t i=0; i<digMult; i++){
677 // Float_t digitAmp = cells.GetCellAmplitude(digID[i]) ;
678 fGeom->RelPosCellInSModule(digID[i], xd, yd, zd);
680 fGeom->GetCellIndex(digID[i],iSM,iT,iIp,iIe);
681 //Gives SuperModule and Tower numbers
684 //*********************************************
686 // fSM[iSupMod]->RegisterCluster(iSM,energy,x,y,z);
688 fSM[iSupMod]->RegisterCluster(iSM,eneInt,x,y,z);
690 } // end cluster loop
693 //______________________________________________________________________________
694 AliEveEMCALSModuleData* AliEveEMCALData::GetSModuleData(Int_t sm)
697 // Return super module data
700 if (sm < 0 || sm > fNsm)
702 printf("The number of super modules must be lower or equal to %d",fNsm);
709 //______________________________________________________________________________
710 void AliEveEMCALData::LoadRaw()
713 // Get raw information
716 // To be implemented !