[u/mrichter/AliRoot.git] / EMCAL / macros / TestEMCALData.C

//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);

  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);
   TTree *sTree = rl->GetTreeS("EMCAL",false);
   TTree *cTree = rl->GetTreeR("EMCAL",false);
   
   TObjArray *harr=NULL;
   TBranch *hbranch=hTree->GetBranch("EMCAL");
   hbranch->SetAddress(&harr);
   
   TObjArray *darr=NULL;
   TBranch *dbranch=dTree->GetBranch("EMCAL");
   dbranch->SetAddress(&darr);
   
   TObjArray *sarr=NULL;
   TBranch *sbranch=sTree->GetBranch("EMCAL");
   sbranch->SetAddress(&sarr);
   
   TObjArray *carr=NULL;
   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){
	 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){
	 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){
	 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){
	 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();
  
}
Commit	Line	Data
66c3bd42	1	//Do a fast test of the different EMCAL data objects
405965e0	2	void TestEMCALData() {
66c3bd42	3
405965e0	4	TH1F* hE = new TH1F("hEmcalHits", "Hits energy distribution in EMCAL", 100, 0., 10.) ;
	5	hE->Sumw2() ;
	6	TH1I* hM = new TH1I("hEmcalHitsMul", "Hits multiplicity distribution in EMCAL", 500, 0., 10000) ;
	7	hM->Sumw2() ;
66c3bd42	8
405965e0	9	TH1I * dA = new TH1I("hEmcalDigits", "Digits amplitude distribution in EMCAL", 500, 0, 5000) ;
	10	dA->Sumw2() ;
	11	TH1I * dM = new TH1I("hEmcalDigitsMul", "Digits multiplicity distribution in EMCAL", 500, 0, 1000) ;
	12	dM->Sumw2() ;
66c3bd42	13
405965e0	14	TH1F * sE = new TH1F("hEmcalSDigits", "SDigits energy distribution in EMCAL", 200, 0, 1000) ;
	15	sE->Sumw2() ;
	16	TH1I * sM = new TH1I("hEmcalSDigitsMul", "SDigits multiplicity distribution in EMCAL", 500, 0, 1000) ;
	17	sM->Sumw2() ;
66c3bd42	18
405965e0	19	TH1F * cE = new TH1F("hEmcalRecPoints", "RecPoints energy distribution in EMCAL", 100, 0, 20) ;
	20	cE->Sumw2() ;
	21	TH1I * cM = new TH1I("hEmcalRecPointsMul", "RecPoints multiplicity distribution in EMCAL", 500, 0, 1000) ;
	22	cM->Sumw2() ;
66c3bd42	23
405965e0	24	AliRunLoader *rl = AliRunLoader::Open("galice.root",AliConfig::GetDefaultEventFolderName(),"read");
	25	if (rl == 0x0)
	26	cout<<"Can not instatiate the Run Loader"<<endl;
66c3bd42	27
405965e0	28	AliEMCALLoader emcalLoader = dynamic_cast<AliEMCALLoader>
405965e0	29	(rl->GetDetectorLoader("EMCAL"));
66c3bd42	30
	31	//Load Hits
	32	rl->LoadHits("EMCAL");
	33	//Load Digits
	34	rl->LoadDigits("EMCAL");
	35	//Load SDigits
	36	rl->LoadSDigits("EMCAL");
	37	//Load RecPoints
	38	rl->LoadRecPoints("EMCAL");
	39
	40
	41	//one way to do it that works for all three
	42	AliEMCALHit* hit;
	43	AliEMCALDigit* dig;
	44	AliEMCALRecPoint* rp;
	45
	46	rl->GetEvent(0);
	47
	48	TClonesArray *hits = 0;
	49	//Fill array of digits
	50	TClonesArray *sdigits = emcalLoader->SDigits();
	51	//Fill array of digits
	52	TClonesArray *digits = emcalLoader->Digits();
	53	//Fill array of clusters
	54	//TObjArray *clusters = emcalLoader->RecPoints(); //It should work, need to FIX
	55	TObjArray *clusters = 0;
	56
	57	//Get hits from the list
	58
	59	//Hits are stored in different branches in the hits Tree,
	60	//first get the branch and then access the hits in the branch
	61	Int_t nHit = 0;
	62	TTree *treeH = emcalLoader->TreeH();
	63	if (treeH) {
	64	// TreeH exists, get the branch
	65	Int_t nTrack = treeH->GetEntries(); // TreeH has array of hits for every primary
	66	TBranch * branchH = treeH->GetBranch("EMCAL");
	67	branchH->SetAddress(&hits);
	68	//Now get the hits in this branch
	69	for (Int_t iTrack = 0; iTrack < nTrack; iTrack++) {
	70	branchH->GetEntry(iTrack);
	71	Int_t nHit = hits->GetEntriesFast();
	72	for(Int_t ihit = 0; ihit< nHit;ihit++){
	73	hit = static_cast<AliEMCALHit *>hits->At(ihit);
	74	if(hit != 0){
	75	nHit++;
	76	hE->Fill(hit->GetEnergy());
	77	cout<<"Hit Info "<<hit->GetId()<<" ELoss "<<hit->GetEnergy()<<endl;
	78	}//hit?
	79	}//hit loop
	80	}// track loop
	81	}//treeH?
	82
	83	hM->Fill(nHit);
	84
	85	//Get digits from the list
	86	if(sdigits){
	87	sM->Fill(sdigits->GetEntries());
	88	for(Int_t isdig = 0; isdig< sdigits->GetEntries();isdig++){
	89	//cout<<">> idig "<<idig<<endl;
	90	dig = static_cast<AliEMCALDigit *>(sdigits->At(isdig)) ;
	91
	92	if(dig != 0){
	93	sE->Fill(dig->GetAmplitude());
94	cout << "SDigit info " << dig->GetId() << " " << dig->GetAmplitude() << endl;
95	}
96	}
97	}
98	else printf("No Sdigits available\n");
99
100
101	//Get digits from the list
102	if(digits){
103	dM->Fill(digits->GetEntries());
104	for(Int_t idig = 0; idig< digits->GetEntries();idig++){
105	//cout<<">> idig "<<idig<<endl;
106	dig = static_cast<AliEMCALDigit *>(digits->At(idig)) ;
107
108	if(dig != 0){
109	dA->Fill(dig->GetAmplitude());
110	cout << "Digit info " << dig->GetId() << " " << dig->GetAmplitude() << endl;
111	}
112	}
113	}
114	else printf("No digits available\n");
115
116	//Get clusters from the list
117	TTree *treeR = emcalLoader->TreeR();
118	TBranch * branchR = treeR->GetBranch("EMCALECARP");
119	branchR->SetAddress(&clusters);
120	branchR->GetEntry(0);
121
122	if(clusters){
123	cM->Fill(clusters->GetEntries());
124	for(Int_t iclu = 0; iclu< clusters->GetEntries();iclu++){
125	//cout<<">> idig "<<idig<<endl;
126	rp = static_cast<AliEMCALRecPoint *>(clusters->At(iclu)) ;
127
128	if(rp != 0){
129	cE->Fill(rp->GetEnergy());
130	cout << "RecPoint info " << rp->GetAbsId(0) << " " << rp->GetEnergy() << endl;
131	}
132	}
133	}else printf("No recpoints available\n");
134
135	/*
405965e0	136	//another way to do it
	137	TTree *hTree = rl->GetTreeH("EMCAL",false);
	138	TTree *dTree = rl->GetTreeD("EMCAL",false);
	139	TTree *sTree = rl->GetTreeS("EMCAL",false);
	140	TTree *cTree = rl->GetTreeR("EMCAL",false);
	141
	142	TObjArray *harr=NULL;
	143	TBranch *hbranch=hTree->GetBranch("EMCAL");
	144	hbranch->SetAddress(&harr);
	145
	146	TObjArray *darr=NULL;
	147	TBranch *dbranch=dTree->GetBranch("EMCAL");
	148	dbranch->SetAddress(&darr);
	149
	150	TObjArray *sarr=NULL;
	151	TBranch *sbranch=sTree->GetBranch("EMCAL");
	152	sbranch->SetAddress(&sarr);
	153
	154	TObjArray *carr=NULL;
	155	TBranch *cbranch=cTree->GetBranch("EMCALECARP");
	156	cbranch->SetAddress(&carr);
	157
66c3bd42	158
405965e0	159	if(hbranch->GetEvent(0)) {
66c3bd42	160	for(Int_t ih = 0; ih < harr->GetEntriesFast(); ih++) {
	161	hM->Fill(harr->GetEntriesFast());
	162	AliEMCALHit* hit =(AliEMCALHit*)harr->UncheckedAt(ih);
	163	if(hit != 0){
405965e0	164	hE->Fill(hit->GetEnergy());
405965e0	165	cout << "Hit info " << hit->GetId() << " " << hit->GetEnergy()*10.5 << endl;
405965e0	166	}
66c3bd42	167	}
	168	}
	169
405965e0	170	if(dbranch->GetEvent(0)) {
66c3bd42	171	for(Int_t id = 0; id < darr->GetEntriesFast(); id++) {
	172	dM->Fill(darr->GetEntriesFast());
	173	AliEMCALDigit* dig =(AliEMCALDigit*)darr->UncheckedAt(id);
	174	if(dig != 0){
405965e0	175	dA->Fill(dig->GetAmp());
405965e0	176	cout << "Digit info " << dig->GetId() << " " << dig->GetAmp() << endl;
405965e0	177	}
66c3bd42	178	}
	179	}
	180
405965e0	181	if(sbranch->GetEvent(0)) {
66c3bd42	182	for(Int_t id = 0; id < sarr->GetEntriesFast(); id++) {
	183	sM->Fill(sarr->GetEntriesFast());
	184	AliEMCALDigit* sdig =(AliEMCALDigit*)sarr->UncheckedAt(id);
	185	if(sdig != 0){
405965e0	186	sE->Fill(sdig->GetAmp()/1.e+6);
405965e0	187	cout << "SDigit info " << sdig->GetId() << " " << sdig->GetAmp()/1.e+6 << endl;
405965e0	188	}
66c3bd42	189	}
	190	}
	191
405965e0	192	if(cbranch->GetEvent(0)) {
66c3bd42	193	for(Int_t ic = 0; ic < carr->GetEntriesFast(); ic++) {
	194	cE->Fill(carr->GetEntriesFast());
	195	AliEMCALRecPoint* rp =(AliEMCALRecPoint*)carr->UncheckedAt(ic);
	196	if(rp != 0){
405965e0	197	cE->Fill(rp->GetEnergy());
405965e0	198	cout << "RecPoint info " << rp->GetAbsId() << " " << rp->GetEnergy() << endl;
66c3bd42	199	}
66c3bd42	200	}
405965e0	201	}
	202
	203	*/
66c3bd42	204
	205	TCanvas *chits = new TCanvas("chits","Hits",20,20,800,400);
	206	chits->Divide(2,1);
	207	chits->cd(1);
	208	hE->Draw();
	209	chits->cd(2);
	210	hM->Draw();
	211
	212	TCanvas *cdig = new TCanvas("cdig","Digits",20,40,800,400);
	213	cdig->Divide(2,1);
	214	cdig->cd(1);
	215	dA->Draw();
	216	cdig->cd(2);
	217	dM->Draw();
	218
	219
	220	TCanvas *csdig = new TCanvas("csdig","SDigits",20,60,800,400);
	221	csdig->Divide(2,1);
	222	csdig->cd(1);
	223	sE->Draw();
	224	csdig->cd(2);
	225	sM->Draw();
	226
	227	TCanvas *cclu = new TCanvas("cclu","Clusters",20,60,800,400);
	228	cclu->Divide(2,1);
	229	cclu->cd(1);
	230	cE->Draw();
	231	cclu->cd(2);
	232	cM->Draw();
	233
405965e0	234	}