update example macros
authorgconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 20 Dec 2010 11:35:33 +0000 (11:35 +0000)
committergconesab <gconesab@f7af4fe6-9843-0410-8265-dc069ae4e863>
Mon, 20 Dec 2010 11:35:33 +0000 (11:35 +0000)
EMCAL/macros/TestEMCALData.C
EMCAL/macros/TestEMCALHit.C [new file with mode: 0644]
EMCAL/macros/TestEMCALRecPoint.C [new file with mode: 0755]

index d85a9c7..8a34907 100644 (file)
-
-
-
+//Do a fast test of the different EMCAL data objects
 void TestEMCALData() {
-
+  
   TH1F* hE = new TH1F("hEmcalHits",    "Hits energy distribution in EMCAL",       100, 0., 10.) ;
   hE->Sumw2() ;
   TH1I* hM = new TH1I("hEmcalHitsMul", "Hits multiplicity distribution in EMCAL", 500, 0., 10000) ;
   hM->Sumw2() ;
-
+  
   TH1I * dA = new TH1I("hEmcalDigits",    "Digits amplitude distribution in EMCAL",    500, 0, 5000) ;
   dA->Sumw2() ;
   TH1I * dM = new TH1I("hEmcalDigitsMul", "Digits multiplicity distribution in EMCAL", 500, 0, 1000) ;
   dM->Sumw2() ;
-
+  
   TH1F * sE = new TH1F("hEmcalSDigits",    "SDigits energy distribution in EMCAL",    200, 0, 1000) ;
   sE->Sumw2() ;
   TH1I * sM = new TH1I("hEmcalSDigitsMul", "SDigits multiplicity distribution in EMCAL", 500, 0, 1000) ;
   sM->Sumw2() ;
-
+  
   TH1F * cE = new TH1F("hEmcalRecPoints",    "RecPoints energy distribution in EMCAL",    100, 0, 20) ;
   cE->Sumw2() ;
   TH1I * cM = new TH1I("hEmcalRecPointsMul", "RecPoints multiplicity distribution in EMCAL", 500, 0, 1000) ;
   cM->Sumw2() ;
-
+  
   AliRunLoader *rl = AliRunLoader::Open("galice.root",AliConfig::GetDefaultEventFolderName(),"read");
   if (rl == 0x0)
     cout<<"Can not instatiate the Run Loader"<<endl;
-
+  
   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>
     (rl->GetDetectorLoader("EMCAL"));
-
-   //Load Hits
-   rl->LoadHits("EMCAL");
-   //Load Digits
-   rl->LoadDigits("EMCAL");
-   //Load SDigits
-   rl->LoadSDigits("EMCAL");
-   //Load RecPoints
-   rl->LoadRecPoints("EMCAL");
-
-   
-   //one way to do it that works for all three
-   AliEMCALHit* hit;
-   AliEMCALDigit* dig;
-   AliEMCALRecPoint* rp;
-
-   rl->GetEvent(0);
-   
-   //Fill array of hits                                                                        
-   TClonesArray *hits = emcalLoader->Hits();
-   //Fill array of digits                                                                        
-   TClonesArray *digits = emcalLoader->Digits();
-   //Fill array of clusters                                                                        
-   TClonesArray *clusters = emcalLoader->RecPoints();
-
-   //Get hits from the list                                                                    
-   hM->Fill(hits->GetEntries());
-   for(Int_t ihit = 0; ihit< hits->GetEntries();ihit++){
-     //cout<<">> idig "<<idig<<endl;                                                             
-     hit = static_cast<AliEMCALHit *>(hits->At(ihit)) ;
-
-     if(hit != 0){
-       hE->Fill(hit->GetEnergy());
-       cout << "Hit info " << hit->GetId() << " " << hit->GetEnergy() << endl;
-     }
-   }
-
-   //Get digits from the list                                                                    
-   dM->Fill(digits->GetEntries());
-   for(Int_t idig = 0; idig< digits->GetEntries();idig++){
-     //cout<<">> idig "<<idig<<endl;                                                             
-     dig = static_cast<AliEMCALDigit *>(digits->At(idig)) ;
-
-     if(dig != 0){
-       dA->Fill(dig->GetAmp());
-       cout << "Digit info " << dig->GetId() << " " << dig->GetAmp() << endl;
-     }
-   }
-
-
-   //Get clusters from the list                                                                    
-   cM->Fill(clusters->GetEntries());
-   for(Int_t iclu = 0; iclu< clusters->GetEntries();iclu++){
-     //cout<<">> idig "<<idig<<endl;                                                             
-     rp = static_cast<AliEMCALRecPoint *>(clusters->At(iclu)) ;
-
-     if(rp != 0){
-       cE->Fill(rp->GetEnergy());
-       cout << "RecPoint info " << rp->GetAbsId() << " " << rp->GetEnergy() << endl;
-     }
-   }
-
-
-   /*
+  
+  //Load Hits
+  rl->LoadHits("EMCAL");
+  //Load Digits
+  rl->LoadDigits("EMCAL");
+  //Load SDigits
+  rl->LoadSDigits("EMCAL");
+  //Load RecPoints
+  rl->LoadRecPoints("EMCAL");
+  
+  
+  //one way to do it that works for all three
+  AliEMCALHit* hit;
+  AliEMCALDigit* dig;
+  AliEMCALRecPoint* rp;
+  
+  rl->GetEvent(0);
+
+  TClonesArray *hits = 0;  
+  //Fill array of digits                                                                        
+  TClonesArray *sdigits = emcalLoader->SDigits();
+  //Fill array of digits                                                                        
+  TClonesArray *digits = emcalLoader->Digits();
+  //Fill array of clusters                                                                        
+  //TObjArray *clusters = emcalLoader->RecPoints(); //It should work, need to FIX
+  TObjArray *clusters = 0;
+  
+  //Get hits from the list  
+  
+  //Hits are stored in different branches in the hits Tree, 
+  //first get the branch and then access the hits in the branch
+  Int_t nHit = 0;
+  TTree *treeH = emcalLoader->TreeH(); 
+  if (treeH) {
+    // TreeH exists, get the branch
+    Int_t nTrack = treeH->GetEntries();  // TreeH has array of hits for every primary
+    TBranch * branchH = treeH->GetBranch("EMCAL");
+    branchH->SetAddress(&hits);
+    //Now get the hits in this branch
+    for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
+      branchH->GetEntry(iTrack);
+      Int_t nHit = hits->GetEntriesFast();
+      for(Int_t ihit = 0; ihit< nHit;ihit++){
+        hit = static_cast<AliEMCALHit *>hits->At(ihit);
+        if(hit != 0){
+          nHit++;
+          hE->Fill(hit->GetEnergy());
+          cout<<"Hit Info "<<hit->GetId()<<" ELoss "<<hit->GetEnergy()<<endl;
+        }//hit?
+      }//hit loop
+    }// track loop
+  }//treeH?
+  
+  hM->Fill(nHit);
+  
+  //Get digits from the list    
+  if(sdigits){
+    sM->Fill(sdigits->GetEntries());
+    for(Int_t isdig = 0; isdig< sdigits->GetEntries();isdig++){
+      //cout<<">> idig "<<idig<<endl;                                                             
+      dig = static_cast<AliEMCALDigit *>(sdigits->At(isdig)) ;
+      
+      if(dig != 0){
+        sE->Fill(dig->GetAmplitude());
+        cout << "SDigit info " << dig->GetId() << " " << dig->GetAmplitude() << endl;
+      }
+    }
+  }
+  else printf("No Sdigits available\n");
+  
+  
+  //Get digits from the list    
+  if(digits){
+    dM->Fill(digits->GetEntries());
+    for(Int_t idig = 0; idig< digits->GetEntries();idig++){
+      //cout<<">> idig "<<idig<<endl;                                                             
+      dig = static_cast<AliEMCALDigit *>(digits->At(idig)) ;
+    
+      if(dig != 0){
+        dA->Fill(dig->GetAmplitude());
+        cout << "Digit info " << dig->GetId() << " " << dig->GetAmplitude() << endl;
+      }
+    }
+  }
+  else printf("No digits available\n");
+  
+  //Get clusters from the list 
+  TTree *treeR = emcalLoader->TreeR();
+  TBranch * branchR = treeR->GetBranch("EMCALECARP");  
+  branchR->SetAddress(&clusters);
+  branchR->GetEntry(0);
+  
+  if(clusters){
+    cM->Fill(clusters->GetEntries());
+    for(Int_t iclu = 0; iclu< clusters->GetEntries();iclu++){
+      //cout<<">> idig "<<idig<<endl;                                                             
+      rp = static_cast<AliEMCALRecPoint *>(clusters->At(iclu)) ;
+      
+      if(rp != 0){
+       cE->Fill(rp->GetEnergy());
+       cout << "RecPoint info " << rp->GetAbsId(0) << " " << rp->GetEnergy() << endl;
+      }
+    }
+  }else printf("No recpoints available\n");
+  
+  /*
    //another way to do it
    TTree *hTree = rl->GetTreeH("EMCAL",false);
    TTree *dTree = rl->GetTreeD("EMCAL",false);
@@ -115,81 +155,80 @@ void TestEMCALData() {
    TBranch *cbranch=cTree->GetBranch("EMCALECARP");
    cbranch->SetAddress(&carr);
    
-
+   
    if(hbranch->GetEvent(0)) {
-     for(Int_t ih = 0; ih < harr->GetEntriesFast(); ih++) {
-       hM->Fill(harr->GetEntriesFast());
-       AliEMCALHit* hit =(AliEMCALHit*)harr->UncheckedAt(ih);
-       if(hit != 0){
+   for(Int_t ih = 0; ih < harr->GetEntriesFast(); ih++) {
+   hM->Fill(harr->GetEntriesFast());
+   AliEMCALHit* hit =(AliEMCALHit*)harr->UncheckedAt(ih);
+   if(hit != 0){
         hE->Fill(hit->GetEnergy());
         cout << "Hit info " << hit->GetId() << " " << hit->GetEnergy()*10.5 << endl;
-       }
-     }
    }
-
+   }
+   }
+   
    if(dbranch->GetEvent(0)) {
-     for(Int_t id = 0; id < darr->GetEntriesFast(); id++) {
-       dM->Fill(darr->GetEntriesFast());
-       AliEMCALDigit* dig =(AliEMCALDigit*)darr->UncheckedAt(id);
-       if(dig != 0){
+   for(Int_t id = 0; id < darr->GetEntriesFast(); id++) {
+   dM->Fill(darr->GetEntriesFast());
+   AliEMCALDigit* dig =(AliEMCALDigit*)darr->UncheckedAt(id);
+   if(dig != 0){
         dA->Fill(dig->GetAmp());
         cout << "Digit info " << dig->GetId() << " " << dig->GetAmp() << endl;
-       }
-     }
    }
-
+   }
+   }
+   
    if(sbranch->GetEvent(0)) {
-     for(Int_t id = 0; id < sarr->GetEntriesFast(); id++) {
-       sM->Fill(sarr->GetEntriesFast());
-       AliEMCALDigit* sdig =(AliEMCALDigit*)sarr->UncheckedAt(id);
-       if(sdig != 0){
+   for(Int_t id = 0; id < sarr->GetEntriesFast(); id++) {
+   sM->Fill(sarr->GetEntriesFast());
+   AliEMCALDigit* sdig =(AliEMCALDigit*)sarr->UncheckedAt(id);
+   if(sdig != 0){
         sE->Fill(sdig->GetAmp()/1.e+6);
         cout << "SDigit info " << sdig->GetId() << " " << sdig->GetAmp()/1.e+6 << endl;
-       }
-     }
    }
-
+   }
+   }
+   
    if(cbranch->GetEvent(0)) {
-     for(Int_t ic = 0; ic < carr->GetEntriesFast(); ic++) {
-       cE->Fill(carr->GetEntriesFast());
-       AliEMCALRecPoint* rp =(AliEMCALRecPoint*)carr->UncheckedAt(ic);
-       if(rp != 0){
+   for(Int_t ic = 0; ic < carr->GetEntriesFast(); ic++) {
+   cE->Fill(carr->GetEntriesFast());
+   AliEMCALRecPoint* rp =(AliEMCALRecPoint*)carr->UncheckedAt(ic);
+   if(rp != 0){
         cE->Fill(rp->GetEnergy());
         cout << "RecPoint info " << rp->GetAbsId() << " " << rp->GetEnergy() << endl;
-       }
-     }
+   }
+   }
    }
    
    */
-       
-   TCanvas *chits = new TCanvas("chits","Hits",20,20,800,400);
-   chits->Divide(2,1);
-   chits->cd(1);
-   hE->Draw();
-   chits->cd(2);
-   hM->Draw();
-       
-   TCanvas *cdig = new TCanvas("cdig","Digits",20,40,800,400);
-   cdig->Divide(2,1);
-   cdig->cd(1);
-   dA->Draw();
-   cdig->cd(2);
-   dM->Draw();
-
-   /*
-   TCanvas *csdig = new TCanvas("csdig","SDigits",20,60,800,400);
-   csdig->Divide(2,1);
-   csdig->cd(1);
-   sE->Draw();
-   csdig->cd(2);
-   sM->Draw();
-   */
-
-   TCanvas *cclu = new TCanvas("cclu","Clusters",20,60,800,400);
-   cclu->Divide(2,1);
-   cclu->cd(1);
-   cE->Draw();
-   cclu->cd(2);
-   cM->Draw();
-
+  
+  TCanvas *chits = new TCanvas("chits","Hits",20,20,800,400);
+  chits->Divide(2,1);
+  chits->cd(1);
+  hE->Draw();
+  chits->cd(2);
+  hM->Draw();
+  
+  TCanvas *cdig = new TCanvas("cdig","Digits",20,40,800,400);
+  cdig->Divide(2,1);
+  cdig->cd(1);
+  dA->Draw();
+  cdig->cd(2);
+  dM->Draw();
+  
+  
+  TCanvas *csdig = new TCanvas("csdig","SDigits",20,60,800,400);
+  csdig->Divide(2,1);
+  csdig->cd(1);
+  sE->Draw();
+  csdig->cd(2);
+  sM->Draw();
+  
+  TCanvas *cclu = new TCanvas("cclu","Clusters",20,60,800,400);
+  cclu->Divide(2,1);
+  cclu->cd(1);
+  cE->Draw();
+  cclu->cd(2);
+  cM->Draw();
+  
 }
diff --git a/EMCAL/macros/TestEMCALHit.C b/EMCAL/macros/TestEMCALHit.C
new file mode 100644 (file)
index 0000000..01585ab
--- /dev/null
@@ -0,0 +1,106 @@
+// Test Macro, shows how to load Hits and Geometry, and how can we get 
+// some of the parameters and variables.
+// Author: Gustavo Conesa
+
+void TestEMCALHit()
+{
+  
+  // Getting EMCAL Detector and Geometry.
+  
+  AliRunLoader *rl = AliRunLoader::Open("galice.root",AliConfig::GetDefaultEventFolderName(),"read");
+  
+  if (rl == 0x0)
+    cout<<"Can not instatiate the Run Loader"<<endl;
+  
+  rl->LoadgAlice();//Needed to get geometry
+  
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>
+    (rl->GetDetectorLoader("EMCAL"));
+  
+  TGeoManager::Import("geometry.root");
+  
+  AliRun * alirun   = rl->GetAliRun(); // Needed to get Geometry
+  AliEMCALGeometry * geom ;
+  if(alirun){
+    AliEMCAL * emcal  = (AliEMCAL*)alirun->GetDetector("EMCAL");
+    geom = emcal->GetGeometry();
+  }
+  
+  if (geom == 0) cout<<"Did not get geometry from EMCALLoader"<<endl;
+  //else   geom->PrintGeometry();
+  
+  //Load Hits
+  rl->LoadHits("EMCAL");
+  
+  //Get maximum number of events
+  Int_t maxevent =  rl->GetNumberOfEvents();
+  cout<<"Number of events "<<maxevent<<endl;
+  //maxevent = 8000 ;
+  
+  
+  AliEMCALHit * hit;
+  TClonesArray *hits = 0;
+  
+  for (Int_t iEvent=0; iEvent<maxevent; iEvent++)
+    {
+      //cout <<  " ======> Event " << iEvent <<endl ;  
+      //Load Event
+      rl->GetEvent(iEvent);
+      Float_t elos=-1;
+      Float_t time  = -1 ;
+      Int_t id      = -1 ;
+      Int_t iSupMod =  0 ;
+      Int_t iTower  =  0 ;
+      Int_t iIphi   =  0 ;
+      Int_t iIeta   =  0 ;
+      Int_t iphi    =  0 ;
+      Int_t ieta    =  0 ;
+      
+      //Fill array of hits
+      cout <<  " ======> Event " << iEvent << endl; 
+      
+      
+      //Get hits from the list      
+      
+      //Hits are stored in different branches in the hits Tree, 
+      //first get the branch and then access the hits in the branch
+      TTree *treeH = emcalLoader->TreeH();     
+      if (treeH) {
+       // TreeH exists, get the branch
+       Int_t nTrack = treeH->GetEntries();  // TreeH has array of hits for every primary
+       TBranch * branchH = treeH->GetBranch("EMCAL");
+       branchH->SetAddress(&hits);
+       for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
+         branchH->GetEntry(iTrack);
+         //Now get the hits in this branch
+         Int_t nHit = hits->GetEntriesFast();
+         for(Int_t ihit = 0; ihit< nHit;ihit++){
+           hit = static_cast<AliEMCALHit *>hits->At(ihit);//(hits->At(ihit)) ;
+           
+           if(hit != 0){
+             id   = hit->GetId() ; //cell (hit) label
+             elos = hit->GetEnergy(); //amplitude in cell (hit)
+             time = hit->GetTime();//time of creation of hit after collision
+             
+             cout<<"Hit ID "<<id<<" ELoss "<<elos;
+             
+             //Geometry methods  
+             if(geom){
+               geom->GetCellIndex(id,iSupMod,iTower,iIphi,iIeta); 
+               //Gives SuperModule and Tower numbers
+               geom->GetCellPhiEtaIndexInSModule(iSupMod,iTower,
+                                                 iIphi, iIeta,iphi,ieta);
+               //Gives label of cell in eta-phi position per each supermodule
+               // cout<< "SModule "<<iSupMod<<"; Tower "<<iTower <<"; Eta "<<iIeta
+               //<<"; Phi "<<iIphi<<"; Cell Eta "<<ieta<<"; Cell Phi "<<iphi<<endl;
+               cout<< ";  SModule "<<iSupMod<<"; Cell Eta "<<ieta<<"; Cell Phi "<<iphi<<endl;
+             }//geom?
+           }//hit?
+           else
+             cout<<"Hit pointer 0x0"<<endl;
+         }//hit loop
+       }// track loop
+      }//treeH?
+    }//event loop
+}
+
diff --git a/EMCAL/macros/TestEMCALRecPoint.C b/EMCAL/macros/TestEMCALRecPoint.C
new file mode 100755 (executable)
index 0000000..1b630da
--- /dev/null
@@ -0,0 +1,95 @@
+// Test Macro, shows how to load RecPoints, and how can we get 
+// some of the parameters and variables.
+// Author: Gustavo Conesa
+
+void TestEMCALRecPoint()
+{
+  
+  // Getting EMCAL Detector and Geometry.
+  
+  AliRunLoader *rl = AliRunLoader::Open("galice.root",AliConfig::GetDefaultEventFolderName(),"read");
+  
+  if (rl == 0x0)
+    cout<<"Can not instatiate the Run Loader"<<endl;
+  
+  rl->LoadgAlice();//Needed to get geometry
+  
+  AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>
+    (rl->GetDetectorLoader("EMCAL"));
+  
+  //TGeoManager::Import("geometry.root");
+  AliGeomManager::LoadGeometry("geometry.root");  
+  AliRun * alirun   = rl->GetAliRun(); // Needed to get Geometry
+  AliEMCALGeometry * geom ;
+  if(alirun){
+    AliEMCAL * emcal  = (AliEMCAL*)alirun->GetDetector("EMCAL");
+    geom = emcal->GetGeometry();
+  }
+  else {
+    cout<<"alirun not available, instantiate"<<endl;
+    geom =  AliEMCALGeometry::GetInstance("EMCAL_COMPLETE") ;  
+  } 
+  
+  //Load RecPoints
+  rl->LoadRecPoints("EMCAL");
+  //Get maximum number of events
+  Int_t maxevent =  rl->GetNumberOfEvents();
+  cout<<"Number of events "<<maxevent<<endl;
+  
+  Int_t iEvent     = -1 ;
+  Int_t iprim      = -1 ;
+  Float_t energy   = -1 ;
+  TVector3 gpos ;
+  
+  AliEMCALRecPoint * rp;
+  for ( iEvent=0; iEvent<maxevent; iEvent++)
+  {
+    cout <<  " ======> Event " << iEvent << endl ;
+    //Load Event
+    rl->GetEvent(iEvent);
+    //      AliStack *sta=rl->Stack();
+    //Fill array of digits
+    TObjArray *rpoints ;//= emcalLoader->RecPoints();
+    
+    TTree *treeR = emcalLoader->TreeR();
+    TBranch * branchR = treeR->GetBranch("EMCALECARP");        
+    branchR->SetAddress(&rpoints);
+    branchR->GetEntry(0);
+    
+    if(!rpoints->GetEntries()) continue;
+    cout <<  " ======> Event " << iEvent << endl ;
+    
+    cout<<">> Entries "<<rpoints->GetEntries()<<endl;
+    
+    //Get recpoints  from the list      
+    for(Int_t irp = 0; irp< rpoints->GetEntries();irp++){
+      rp = static_cast<AliEMCALRecPoint *>(rpoints->At(irp)) ;
+      //rp = emcalLoader->RecPoint(irp);
+      if(rp != 0){
+        {
+          energy  = rp->GetEnergy(); //cluster energy
+          cout<<"Energy "<<energy<<" PointEnergy "<<rp->GetPointEnergy()<<endl;
+          rp->GetGlobalPosition(gpos);//global ALICE xyz position
+          Int_t  primMult  = 0;
+          //Int_t *primInts =  rp->GetPrimaries(primMult);
+          //for (Int_t ipr=0; ipr<primMult; ipr++) 
+          //   cout<<"primlist "<<ipr<<" index "<< primInts[ipr]<<endl;
+          //iprim=rp->GetPrimaryIndex() ;
+          //    cout<<"Selected primary index "<<iprim<<endl;
+          //  if(iprim>0){
+          //TParticle *primary=sta->Particle(iprim);
+          //cout<<"Primary phi "<<primary->Phi()*180/TMath::Pi()<<" Reconstructed phi "<<gpos.Phi()*180/TMath::Pi()<<"; Energy "<<primary->Energy()<<endl;
+          //h->Fill((gpos.Phi()-primary->Phi())*180/TMath::Pi());
+          cout<<"rec point "<<irp<<"; Energy "<<energy<<" Phi "<<gpos.Phi()*180/TMath::Pi()<<" Eta "<<gpos.Eta()<<" iprim "<<iprim<<endl; 
+          //}
+          Float_t lmb[2];
+          rp->GetElipsAxis(lmb);
+          cout<<"lmb0 "<<lmb[0]<<" lmb1 "<<lmb[1]<<endl;
+        }
+      }
+      else
+        cout<<"recpoint pointer 0x0"<<endl;
+    } 
+  }
+}
+