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"
22 #include "AliEMCALLoader.h"
23 #include "AliESDVertex.h"
24 #include "AliEMCALHit.h"
25 #include "AliEMCALDigit.h"
27 #include "AliEMCALRecPoint.h"
28 #include "AliESDCaloCells.h"
29 #include "AliESDCaloCluster.h"
31 #include "AliEveEMCALData.h"
32 #include "AliEveEMCALSModuleData.h"
41 class AliEMCALGeometry;
42 class AliEveEMCALSModule;
45 ClassImp(AliEveEMCALData)
47 //______________________________________________________________________________
48 AliEveEMCALData::AliEveEMCALData():
71 // for(Int_t i=0; i<12; i++)
73 // for(Int_t i=0; i<10; i++)
75 // fSMhalf[0] = fSMhalf[1] = 0x0;
82 //______________________________________________________________________________
83 AliEveEMCALData::AliEveEMCALData(AliRunLoader* rl, TGeoNode* node, TGeoHMatrix* m):
106 // for(Int_t i=0; i<12; i++)
108 // for(Int_t i=0; i<10; i++)
110 // fSMhalf[0] = fSMhalf[1] = 0x0;
117 //______________________________________________________________________________
118 AliEveEMCALData::~AliEveEMCALData()
124 DeleteSuperModules();
126 // delete fEmcal; // deleted by run-loader
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* const tree)
179 // Set digit-tree to be used for digit retrieval.
180 // Data is loaded on demand.
187 //______________________________________________________________________________
188 void AliEveEMCALData::SetESD(AliESDEvent* const esd)
197 //______________________________________________________________________________
198 void AliEveEMCALData::SetNode(TGeoNode* const node)
207 //______________________________________________________________________________
208 void AliEveEMCALData::InitEMCALGeom(AliRunLoader* const 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* const 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* const emcl)
356 // Get hit information from EMCAL Loader
362 TClonesArray *hits = 0;//(TClonesArray*)emcl->Hits();
363 TTree *treeH = emcl->TreeH();
365 Int_t nTrack = treeH->GetEntries(); // TreeH has array of hits for every primary
366 TBranch * branchH = treeH->GetBranch("EMCAL");
367 //if(fHits)fHits->Clear();
368 branchH->SetAddress(&hits);
369 for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
370 branchH->GetEntry(iTrack);
372 //Get hits from the list
373 for(Int_t ihit = 0; ihit< hits->GetEntries();ihit++){
375 hit = static_cast<AliEMCALHit *>(hits->At(ihit)) ;
378 if(fDebug>1) cout << "Hit info " << hit->GetId() << " " << hit->GetEnergy() << endl;
380 Int_t id = hit->GetId();
381 // These are local coordinates
382 Double_t xl = 0.; Double_t yl = 0.; Double_t zl = 0.;
383 // Get global coordinates
384 Double_t x = hit->X();
385 Double_t y = hit->Y();
386 Double_t z = hit->Z();
387 Double_t amp = hit->GetEnergy();
390 GetGeomInfo(id,iSupMod,xl,yl,zl);
391 fSM[iSupMod]->RegisterHit(id,iSupMod,amp,x,y,z);
399 //______________________________________________________________________________
400 void AliEveEMCALData::LoadDigits(TTree *t)
403 // Get digit information from RunLoader
406 TClonesArray *digits = 0;
407 t->SetBranchAddress("EMCAL", &digits);
410 Int_t nEnt = digits->GetEntriesFast();
411 cout << "nEnt: " << nEnt << endl;
414 // Double_t amp = -1 ;
415 Double_t ampFlo = -1 ;
420 for (Int_t idig = 0; idig < nEnt; idig++)
422 dig = static_cast<AliEMCALDigit *>(digits->At(idig));
425 id = dig->GetId() ; //cell (digit) label
427 ampFlo = dig->GetAmplitude(); //amplitude in cell (digit)
429 // amp = ampFlo*0.0153; // To be modified with correct OCDB conversion
431 GetGeomInfo(id,iSupMod,x,y,z);
434 // fSM[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
435 // // fSM[iSupMod]->SaveDigit(dig);
436 // // if(iSupMod<fNsmfull) fSMfull[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
437 // // if(iSupMod>fNsmfull) fSMhalf[iSupMod-10]->RegisterDigit(id,iSupMod,amp,x,y,z);
438 fSM[iSupMod]->RegisterDigit(id,iSupMod,ampFlo,x,y,z);
439 // fSM[iSupMod]->SaveDigit(dig);
440 // if(iSupMod<fNsmfull) fSMfull[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
441 // if(iSupMod>fNsmfull) fSMhalf[iSupMod-10]->RegisterDigit(id,iSupMod,amp,x,y,z);
444 cout << "Digit object empty" << endl;
448 cout << "after loop on digits !" << endl;
451 //______________________________________________________________________________
452 void AliEveEMCALData::LoadDigitsFromEMCALLoader(AliEMCALLoader* const emcl)
456 // Get digit information from EMCAL Loader
461 //Fill array of digits
462 TClonesArray *digits = (TClonesArray*)emcl->Digits();
464 //Get digits from the list
466 // Double_t amp = -1 ;
467 Double_t ampFlo = -1 ;
472 for(Int_t idig = 0; idig< digits->GetEntries();idig++){
474 dig = static_cast<AliEMCALDigit *>(digits->At(idig)) ;
477 if(fDebug>1) cout << "Digit info " << dig->GetId() << " " << dig->GetAmplitude() << endl;
478 id = dig->GetId() ; //cell (digit) label
480 ampFlo = dig->GetAmplitude(); //amplitude in cell (digit)
482 // amp = ampFlo*0.0153.; // To be modified with correct OCDB conversion
484 GetGeomInfo(id,iSupMod,x,y,z);
487 // fSM[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
489 fSM[iSupMod]->RegisterDigit(id,iSupMod,ampFlo,x,y,z);
492 cout << "Digit object empty" << endl;
495 } // end loop on digits
499 //______________________________________________________________________________
500 void AliEveEMCALData::LoadDigitsFromESD()
503 // Get digit information from esd
506 AliESDCaloCells &cells= *(fESD->GetEMCALCells());
507 Int_t ncell = cells.GetNumberOfCells() ;
511 // Extract digit information from the ESDs
512 for (Int_t icell= 0; icell < ncell; icell++)
514 Int_t id = cells.GetCellNumber(icell);
516 Double_t ampFlo = cells.GetAmplitude(icell);
518 // Double_t amp = ampFlo*0.0153; // To be modified with correct OCDB conversion
520 GetGeomInfo(id,iSupMod,x,y,z);
523 // fSM[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
524 // if(iSupMod<fNsmfull) fSMfull[iSupMod]->RegisterDigit(id,iSupMod,amp,x,y,z);
525 // if(iSupMod>fNsmfull) fSMhalf[iSupMod-10]->RegisterDigit(id,iSupMod,amp,x,y,z);
527 fSM[iSupMod]->RegisterDigit(id,iSupMod,ampFlo,x,y,z);
528 if(iSupMod<fNsmfull) fSMfull[iSupMod]->RegisterDigit(id,iSupMod,ampFlo,x,y,z);
529 if(iSupMod>fNsmfull) fSMhalf[iSupMod-10]->RegisterDigit(id,iSupMod,ampFlo,x,y,z);
534 //______________________________________________________________________________
535 void AliEveEMCALData::LoadRecPoints(TTree* const t)
538 // Get rec point information from RunLoader
541 //*************************************************
542 // To be improved !!!!!
543 // Size and shape of cluster to be implemented
545 //*************************************************
548 TObjArray *carr=NULL;
549 TBranch *cbranch=t->GetBranch("EMCALECARP");
550 cbranch->SetAddress(&carr);
552 if(cbranch->GetEvent(0)) {
553 for(Int_t ic = 0; ic < carr->GetEntriesFast(); ic++) {
554 AliEMCALRecPoint* rp =(AliEMCALRecPoint*)carr->UncheckedAt(ic);
556 if(fDebug>1) cout << "RecPoint info " << rp->GetAbsId() << " " << rp->GetEnergy() << endl;
557 Int_t iSupMod = rp->GetSuperModuleNumber();
559 Double_t amp = (Double_t)rp->GetEnergy();
561 Double_t ampFlo = amp/0.0153; // To be modified with correct OCDB conversion
563 rp->GetLocalPosition(lpos);
566 // fSM[iSupMod]->RegisterCluster(iSupMod,amp,lpos[0],lpos[1],lpos[2]);
568 fSM[iSupMod]->RegisterCluster(iSupMod,ampFlo,lpos[0],lpos[1],lpos[2]);
575 //______________________________________________________________________________
576 void AliEveEMCALData::LoadRecPointsFromEMCALLoader(AliEMCALLoader* const emcl)
579 // Get rec point information from EMCAL Loader
582 //*************************************************
583 // To be improved !!!!!
584 // Size and shape of cluster to be implemented
586 //*************************************************
589 AliEMCALRecPoint* rp;
591 //Fill array of clusters
592 TClonesArray *clusters = (TClonesArray*)emcl->RecPoints();
594 //Get clusters from the list
595 for(Int_t iclu = 0; iclu< clusters->GetEntries();iclu++){
597 rp = static_cast<AliEMCALRecPoint *>(clusters->At(iclu)) ;
600 if(fDebug>1) cout << "RecPoint info " << rp->GetAbsId() << " " << rp->GetEnergy() << endl;
601 Int_t iSupMod = rp->GetSuperModuleNumber();
602 Double_t amp = (Double_t)rp->GetEnergy();
603 Double_t ampFlo = amp/0.0153; // To be modified with correct OCDB conversion
605 rp->GetLocalPosition(lpos);
608 // fSM[iSupMod]->RegisterCluster(iSupMod,amp,lpos[0],lpos[1],lpos[2]);
610 fSM[iSupMod]->RegisterCluster(iSupMod,ampFlo,lpos[0],lpos[1],lpos[2]);
616 //______________________________________________________________________________
617 void AliEveEMCALData::LoadRecPointsFromESD()
620 // Get cluster information from esd
632 // Get reconstructed vertex position
633 AliESDVertex* primVertex =(AliESDVertex*) fESD->GetVertex();
634 Double_t vertexPosition[3] ;
635 primVertex->GetXYZ(vertexPosition) ;
637 //Get the CaloClusters
638 //select EMCAL clusters only
639 TRefArray * caloClusters = new TRefArray();
640 fESD->GetEMCALClusters(caloClusters);
641 Int_t nclus = caloClusters->GetEntries();
642 cout << "nclus: " << nclus << endl;
644 for (Int_t iclus = 0; iclus < nclus; iclus++)
646 AliESDCaloCluster *clus = (AliESDCaloCluster *) caloClusters->At(iclus) ;
647 //Get the cluster info
649 Double_t energy = clus->E() ;
651 // Int_t eneInt = (Int_t)energy*500+0.5;
652 Double_t eneInt = energy/0.0153; // To be modified with correct OCDB conversion
653 Double_t disp = clus->GetDispersion() ;
655 clus->GetPosition(pos) ; // Global position
656 TVector3 vpos(pos[0],pos[1],pos[2]) ;
659 clus->GetMomentum(p4,vertexPosition);
660 p3.SetXYZ(p4[0],p4[1],p4[2]);
661 Double_t eta = p3.Eta();
662 Double_t phi = ( (p3.Phi()) < 0) ? (p3.Phi()) + 2. * TMath::Pi() : (p3.Phi());
664 Int_t mult = clus->GetNCells() ;
666 cout << "In cluster: " << iclus << ", ncells: " << mult << ", energy : " << energy <<
667 ", disp: " << disp << endl;
668 cout << "Cluster " << iclus << ", eta: " << eta << ", phi: " << phi << endl;
672 fGeom->GetAbsCellIdFromEtaPhi(eta,phi,clusId);
674 cout << "Abs Cluster Id: " << clusId << ", xc: " << pos[0] <<
675 ", yc: " << pos[1] << ", zc: " << pos[2] << endl;
678 GetGeomInfo(clusId,iSupMod,x,y,z);
680 //******** Not used yet but will come ********
681 // AliESDCaloCells &cells= *(fESD->GetEMCALCells());
682 Int_t digMult = clus->GetNCells() ;
683 UShort_t *digID = clus->GetCellsAbsId() ;
684 for(Int_t i=0; i<digMult; i++){
685 // Float_t digitAmp = cells.GetCellAmplitude(digID[i]) ;
686 fGeom->RelPosCellInSModule(digID[i], xd, yd, zd);
688 fGeom->GetCellIndex(digID[i],iSM,iT,iIp,iIe);
689 //Gives SuperModule and Tower numbers
692 //*********************************************
694 // fSM[iSupMod]->RegisterCluster(iSM,energy,x,y,z);
696 fSM[iSupMod]->RegisterCluster(iSM,eneInt,x,y,z);
698 } // end cluster loop
701 //______________________________________________________________________________
702 AliEveEMCALSModuleData* AliEveEMCALData::GetSModuleData(Int_t sm)
705 // Return super module data
708 if (sm < 0 || sm > fNsm)
710 printf("The number of super modules must be lower or equal to %d",fNsm);
717 //______________________________________________________________________________
718 void AliEveEMCALData::LoadRaw() const
721 // Get raw information
724 // To be implemented !