2 // Fill containers for visualisation of EMCAL data structures
3 // - read and store MC Hits - read and store digits from esds or runloader
4 // - read and store clusters from esds or runloader
6 // Author: Magali Estienne (magali.estienne@cern.ch)
10 //#include <Riostream.h>
15 #include <TObjArray.h>
16 #include <TRefArray.h>
17 #include <TClonesArray.h>
19 #include <TLorentzVector.h>
21 #include "AliRunLoader.h"
23 #include "AliEMCALLoader.h"
24 #include "AliESDVertex.h"
25 #include "AliEMCALHit.h"
26 #include "AliEMCALDigit.h"
28 #include "AliEMCALRecPoint.h"
29 #include "AliESDCaloCells.h"
30 #include "AliESDCaloCluster.h"
32 #include "AliEveEMCALData.h"
33 #include "AliEveEMCALSModuleData.h"
42 class AliEMCALGeometry;
43 class AliEveEMCALSModule;
47 ClassImp(AliEveEMCALData)
49 //______________________________________________________________________________
50 AliEveEMCALData::AliEveEMCALData():
73 // for(Int_t i=0; i<12; i++)
75 // for(Int_t i=0; i<10; i++)
77 // fSMhalf[0] = fSMhalf[1] = 0x0;
84 //______________________________________________________________________________
85 AliEveEMCALData::AliEveEMCALData(AliRunLoader* rl, TGeoNode* node, TGeoHMatrix* m):
108 // for(Int_t i=0; i<12; i++)
110 // for(Int_t i=0; i<10; i++)
112 // fSMhalf[0] = fSMhalf[1] = 0x0;
119 //______________________________________________________________________________
120 AliEveEMCALData::~AliEveEMCALData()
126 DeleteSuperModules();
128 // delete fEmcal; // deleted by run-loader
135 //______________________________________________________________________________
136 AliEveEMCALData::AliEveEMCALData(const AliEveEMCALData &edata) :
139 fEmcal(edata.fEmcal),
142 fHMatrix(edata.fHMatrix),
146 fNsmfull(edata.fNsmfull),
147 fNsmhalf(edata.fNsmhalf),
149 fSMfull(edata.fSMfull),
150 fSMhalf(edata.fSMhalf),
151 fRunLoader(edata.fRunLoader),
152 fDebug(edata.fDebug),
158 InitEMCALGeom(edata.fRunLoader);
162 //______________________________________________________________________________
163 AliEveEMCALData& AliEveEMCALData::operator=(const AliEveEMCALData &edata)
166 // Assignment operator
169 if (this != &edata) {
177 //______________________________________________________________________________
178 void AliEveEMCALData::SetTree(TTree* const tree)
181 // Set digit-tree to be used for digit retrieval.
182 // Data is loaded on demand.
189 //______________________________________________________________________________
190 void AliEveEMCALData::SetESD(AliESDEvent* const esd)
199 //______________________________________________________________________________
200 void AliEveEMCALData::SetNode(TGeoNode* const node)
209 //______________________________________________________________________________
210 void AliEveEMCALData::InitEMCALGeom(AliRunLoader* const rl)
213 // Set data members for EMCAL geometry
216 fEmcal = (AliEMCAL*) rl->GetAliRun()->GetDetector("EMCAL");
217 fGeom = (AliEMCALGeometry*) fEmcal->GetGeometry();
221 //______________________________________________________________________________
222 void AliEveEMCALData::GetGeomInfo(Int_t id, Int_t &iSupMod, Double_t& x, Double_t& y, Double_t& z)
225 // Get geometrical information from hit/digit/cluster absolute id
233 fGeom->GetCellIndex(id,iSupMod,iTower,iIphi,iIeta);
234 //Gives SuperModule and Tower numbers
235 fGeom->RelPosCellInSModule(id, x, y, z);
239 //______________________________________________________________________________
240 void AliEveEMCALData::CreateAllSModules()
243 // Create all fNsm super modules
246 for(Int_t sm = 0; sm < fNsm; sm++)
251 //______________________________________________________________________________
252 void AliEveEMCALData::CreateSModule(Int_t sm)
255 // Create super-module-data for SM if it does not exist already.
258 if(fSM[sm] == 0) fSM[sm] = new AliEveEMCALSModuleData(sm,fGeom,fNode,fHMatrix);
259 if(fSMfull[sm] == 0 && sm < 10) fSMfull[sm] = new AliEveEMCALSModuleData(sm,fGeom,fNode,fHMatrix);
260 if(fSMhalf[sm-10] == 0 && sm > 10) fSMhalf[sm-10] = new AliEveEMCALSModuleData(sm,fGeom,fNode,fHMatrix);
263 //______________________________________________________________________________
264 void AliEveEMCALData::DropAllSModules()
267 // Drop data of all existing sectors.
270 for (Int_t sm = 0; sm < fNsm; sm++) {
276 //______________________________________________________________________________
277 void AliEveEMCALData::DeleteSuperModules()
280 // Delete all super module data
283 for (Int_t sm = 0; sm < fNsm; sm++)
289 for(Int_t smf = 0; smf < fNsmfull; smf++)
295 for(Int_t smh = 0; smh < fNsmhalf; smh++)
303 //______________________________________________________________________________
304 void AliEveEMCALData::LoadHits(TTree* const t)
307 // Get hit information from RunLoader
311 // These are global coordinates !
313 const char *selection = "";
314 const char *varexp = "fX:fY:fZ";
315 sprintf(form,"EMCAL Hits '%s'", selection);
316 fPoint = new TEvePointSet(form);
318 TEvePointSelector ps(t, fPoint, varexp, selection);
321 if (fPoint->Size() == 0) {
322 Warning("emcal_hits", Form("No hits match '%s'", selection));
328 TObjArray *harr=NULL;
329 TBranch *hbranch=t->GetBranch("EMCAL");
330 hbranch->SetAddress(&harr);
332 if(hbranch->GetEvent(0)) {
333 for(Int_t ih = 0; ih < harr->GetEntriesFast(); ih++) {
334 AliEMCALHit* hit =(AliEMCALHit*)harr->UncheckedAt(ih);
336 if(fDebug>1) cout << "Hit info " << hit->GetId() << " " << hit->GetEnergy() << endl;
337 Int_t id = hit->GetId();
338 // These are local coordinates
339 Double_t xl = 0.; Double_t yl = 0.; Double_t zl = 0.;
340 // Get global coordinates
341 Double_t x = hit->X();
342 Double_t y = hit->Y();
343 Double_t z = hit->Z();
344 Double_t amp = hit->GetEnergy();
347 GetGeomInfo(id,iSupMod,xl,yl,zl);
348 fSM[iSupMod]->RegisterHit(id,iSupMod,amp,x,y,z);
354 //______________________________________________________________________________
355 void AliEveEMCALData::LoadHitsFromEMCALLoader(AliEMCALLoader* const emcl)
358 // Get hit information from EMCAL Loader
364 TClonesArray *hits = 0;//(TClonesArray*)emcl->Hits();
365 TTree *treeH = emcl->TreeH();
367 Int_t nTrack = treeH->GetEntries(); // TreeH has array of hits for every primary
368 TBranch * branchH = treeH->GetBranch("EMCAL");
369 //if(fHits)fHits->Clear();
370 branchH->SetAddress(&hits);
371 for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
372 branchH->GetEntry(iTrack);
374 //Get hits from the list
375 for(Int_t ihit = 0; ihit< hits->GetEntries();ihit++){
377 hit = static_cast<AliEMCALHit *>(hits->At(ihit)) ;
380 if(fDebug>1) cout << "Hit info " << hit->GetId() << " " << hit->GetEnergy() << endl;
382 Int_t id = hit->GetId();
383 // These are local coordinates
384 Double_t xl = 0.; Double_t yl = 0.; Double_t zl = 0.;
385 // Get global coordinates
386 Double_t x = hit->X();
387 Double_t y = hit->Y();
388 Double_t z = hit->Z();
389 Double_t amp = hit->GetEnergy();
392 GetGeomInfo(id,iSupMod,xl,yl,zl);
393 fSM[iSupMod]->RegisterHit(id,iSupMod,amp,x,y,z);
401 //______________________________________________________________________________
402 void AliEveEMCALData::LoadDigits(TTree *t)
405 // Get digit information from RunLoader
408 TClonesArray *digits = 0;
409 t->SetBranchAddress("EMCAL", &digits);
412 Int_t nEnt = digits->GetEntriesFast();
413 cout << "nEnt: " << nEnt << endl;
416 // Double_t amp = -1 ;
417 Double_t ampFlo = -1 ;
422 for (Int_t idig = 0; idig < nEnt; idig++)
424 dig = static_cast<AliEMCALDigit *>(digits->At(idig));
427 id = dig->GetId() ; //cell (digit) label
429 ampFlo = dig->GetAmplitude(); //amplitude in cell (digit)
431 // amp = ampFlo*0.0153; // To be modified with correct OCDB conversion
433 GetGeomInfo(id,iSupMod,x,y,z);
436 // fSM[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
437 // // fSM[iSupMod]->SaveDigit(dig);
438 // // if(iSupMod<fNsmfull) fSMfull[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
439 // // if(iSupMod>fNsmfull) fSMhalf[iSupMod-10]->RegisterDigit(id,iSupMod,amp,x,y,z);
440 fSM[iSupMod]->RegisterDigit(id,iSupMod,ampFlo,x,y,z);
441 // fSM[iSupMod]->SaveDigit(dig);
442 // if(iSupMod<fNsmfull) fSMfull[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
443 // if(iSupMod>fNsmfull) fSMhalf[iSupMod-10]->RegisterDigit(id,iSupMod,amp,x,y,z);
446 cout << "Digit object empty" << endl;
450 cout << "after loop on digits !" << endl;
453 //______________________________________________________________________________
454 void AliEveEMCALData::LoadDigitsFromEMCALLoader(AliEMCALLoader* const emcl)
458 // Get digit information from EMCAL Loader
463 //Fill array of digits
464 TClonesArray *digits = (TClonesArray*)emcl->Digits();
466 //Get digits from the list
468 // Double_t amp = -1 ;
469 Double_t ampFlo = -1 ;
474 for(Int_t idig = 0; idig< digits->GetEntries();idig++){
476 dig = static_cast<AliEMCALDigit *>(digits->At(idig)) ;
479 if(fDebug>1) cout << "Digit info " << dig->GetId() << " " << dig->GetAmplitude() << endl;
480 id = dig->GetId() ; //cell (digit) label
482 ampFlo = dig->GetAmplitude(); //amplitude in cell (digit)
484 // amp = ampFlo*0.0153.; // To be modified with correct OCDB conversion
486 GetGeomInfo(id,iSupMod,x,y,z);
489 // fSM[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
491 fSM[iSupMod]->RegisterDigit(id,iSupMod,ampFlo,x,y,z);
494 cout << "Digit object empty" << endl;
497 } // end loop on digits
501 //______________________________________________________________________________
502 void AliEveEMCALData::LoadDigitsFromESD()
505 // Get digit information from esd
508 AliESDCaloCells &cells= *(fESD->GetEMCALCells());
509 Int_t ncell = cells.GetNumberOfCells() ;
513 // Extract digit information from the ESDs
514 for (Int_t icell= 0; icell < ncell; icell++)
516 Int_t id = cells.GetCellNumber(icell);
518 Double_t ampFlo = cells.GetAmplitude(icell);
520 // Double_t amp = ampFlo*0.0153; // To be modified with correct OCDB conversion
522 GetGeomInfo(id,iSupMod,x,y,z);
525 // fSM[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
526 // if(iSupMod<fNsmfull) fSMfull[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
527 // if(iSupMod>fNsmfull) fSMhalf[iSupMod-10]->RegisterDigit(id,iSupMod,amp,x,y,z);
529 fSM[iSupMod]->RegisterDigit(id,iSupMod,ampFlo,x,y,z);
530 if(iSupMod<fNsmfull) fSMfull[iSupMod]->RegisterDigit(id,iSupMod,ampFlo,x,y,z);
531 if(iSupMod>fNsmfull) fSMhalf[iSupMod-10]->RegisterDigit(id,iSupMod,ampFlo,x,y,z);
536 //______________________________________________________________________________
537 void AliEveEMCALData::LoadRecPoints(TTree* const t)
540 // Get rec point information from RunLoader
543 //*************************************************
544 // To be improved !!!!!
545 // Size and shape of cluster to be implemented
547 //*************************************************
550 TObjArray *carr=NULL;
551 TBranch *cbranch=t->GetBranch("EMCALECARP");
552 cbranch->SetAddress(&carr);
554 if(cbranch->GetEvent(0)) {
555 for(Int_t ic = 0; ic < carr->GetEntriesFast(); ic++) {
556 AliEMCALRecPoint* rp =(AliEMCALRecPoint*)carr->UncheckedAt(ic);
558 if(fDebug>1) cout << "RecPoint info " << rp->GetAbsId() << " " << rp->GetEnergy() << endl;
559 Int_t iSupMod = rp->GetSuperModuleNumber();
561 Double_t amp = (Double_t)rp->GetEnergy();
563 Double_t ampFlo = amp/0.0153; // To be modified with correct OCDB conversion
565 rp->GetLocalPosition(lpos);
568 // fSM[iSupMod]->RegisterCluster(iSupMod,amp,lpos[0],lpos[1],lpos[2]);
570 fSM[iSupMod]->RegisterCluster(iSupMod,ampFlo,lpos[0],lpos[1],lpos[2]);
577 //______________________________________________________________________________
578 void AliEveEMCALData::LoadRecPointsFromEMCALLoader(AliEMCALLoader* const emcl)
581 // Get rec point information from EMCAL Loader
584 //*************************************************
585 // To be improved !!!!!
586 // Size and shape of cluster to be implemented
588 //*************************************************
591 AliEMCALRecPoint* rp;
593 //Fill array of clusters
594 TClonesArray *clusters = (TClonesArray*)emcl->RecPoints();
596 //Get clusters from the list
597 for(Int_t iclu = 0; iclu< clusters->GetEntries();iclu++){
599 rp = static_cast<AliEMCALRecPoint *>(clusters->At(iclu)) ;
602 if(fDebug>1) cout << "RecPoint info " << rp->GetAbsId() << " " << rp->GetEnergy() << endl;
603 Int_t iSupMod = rp->GetSuperModuleNumber();
604 Double_t amp = (Double_t)rp->GetEnergy();
605 Double_t ampFlo = amp/0.0153; // To be modified with correct OCDB conversion
607 rp->GetLocalPosition(lpos);
610 // fSM[iSupMod]->RegisterCluster(iSupMod,amp,lpos[0],lpos[1],lpos[2]);
612 fSM[iSupMod]->RegisterCluster(iSupMod,ampFlo,lpos[0],lpos[1],lpos[2]);
618 //______________________________________________________________________________
619 void AliEveEMCALData::LoadRecPointsFromESD()
622 // Get cluster information from esd
634 // Get reconstructed vertex position
635 AliESDVertex* primVertex =(AliESDVertex*) fESD->GetVertex();
636 Double_t vertexPosition[3] ;
637 primVertex->GetXYZ(vertexPosition) ;
639 //Get the CaloClusters
640 //select EMCAL clusters only
641 TRefArray * caloClusters = new TRefArray();
642 fESD->GetEMCALClusters(caloClusters);
643 Int_t nclus = caloClusters->GetEntries();
644 cout << "nclus: " << nclus << endl;
646 for (Int_t iclus = 0; iclus < nclus; iclus++)
648 AliESDCaloCluster *clus = (AliESDCaloCluster *) caloClusters->At(iclus) ;
649 //Get the cluster info
651 Double_t energy = clus->E() ;
653 // Int_t eneInt = (Int_t)energy*500+0.5;
654 Double_t eneInt = energy/0.0153; // To be modified with correct OCDB conversion
655 Double_t disp = clus->GetDispersion() ;
657 clus->GetPosition(pos) ; // Global position
658 TVector3 vpos(pos[0],pos[1],pos[2]) ;
661 clus->GetMomentum(p4,vertexPosition);
662 p3.SetXYZ(p4[0],p4[1],p4[2]);
663 Double_t eta = p3.Eta();
664 Double_t phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi());
666 Int_t mult = clus->GetNCells() ;
668 cout << "In cluster: " << iclus << ", ncells: " << mult << ", energy : " << energy <<
669 ", disp: " << disp << endl;
670 cout << "Cluster " << iclus << ", eta: " << eta << ", phi: " << phi << endl;
674 fGeom->GetAbsCellIdFromEtaPhi(eta,phi,clusId);
676 cout << "Abs Cluster Id: " << clusId << ", xc: " << pos[0] <<
677 ", yc: " << pos[1] << ", zc: " << pos[2] << endl;
680 GetGeomInfo(clusId,iSupMod,x,y,z);
682 //******** Not used yet but will come ********
683 // AliESDCaloCells &cells= *(fESD->GetEMCALCells());
684 Int_t digMult = clus->GetNCells() ;
685 UShort_t *digID = clus->GetCellsAbsId() ;
686 for(Int_t i=0; i<digMult; i++){
687 // Float_t digitAmp = cells.GetCellAmplitude(digID[i]) ;
688 fGeom->RelPosCellInSModule(digID[i], xd, yd, zd);
690 fGeom->GetCellIndex(digID[i],iSM,iT,iIp,iIe);
691 //Gives SuperModule and Tower numbers
694 //*********************************************
696 // fSM[iSupMod]->RegisterCluster(iSM,energy,x,y,z);
698 fSM[iSupMod]->RegisterCluster(iSM,eneInt,x,y,z);
700 } // end cluster loop
703 //______________________________________________________________________________
704 AliEveEMCALSModuleData* AliEveEMCALData::GetSModuleData(Int_t sm)
707 // Return super module data
710 if (sm < 0 || sm > fNsm)
712 printf("The number of super modules must be lower or equal to %d",fNsm);
719 //______________________________________________________________________________
720 void AliEveEMCALData::LoadRaw() const
723 // Get raw information
726 // To be implemented !