suggestion by Chiara to add TestPreprocessor macro to svn
[u/mrichter/AliRoot.git] / EMCAL / macros / TestEMCALData.C
1 //Do a fast test of the different EMCAL data objects
2 void TestEMCALData() {
3   
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() ;
8   
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() ;
13   
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() ;
18   
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() ;
23   
24   AliRunLoader *rl = AliRunLoader::Open("galice.root",AliConfig::GetDefaultEventFolderName(),"read");
25   if (rl == 0x0)
26     cout<<"Can not instatiate the Run Loader"<<endl;
27   
28   AliEMCALLoader *emcalLoader = dynamic_cast<AliEMCALLoader*>
29     (rl->GetDetectorLoader("EMCAL"));
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   /*
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    
158    
159    if(hbranch->GetEvent(0)) {
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){
164          hE->Fill(hit->GetEnergy());
165          cout << "Hit info " << hit->GetId() << " " << hit->GetEnergy()*10.5 << endl;
166    }
167    }
168    }
169    
170    if(dbranch->GetEvent(0)) {
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){
175          dA->Fill(dig->GetAmp());
176          cout << "Digit info " << dig->GetId() << " " << dig->GetAmp() << endl;
177    }
178    }
179    }
180    
181    if(sbranch->GetEvent(0)) {
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){
186          sE->Fill(sdig->GetAmp()/1.e+6);
187          cout << "SDigit info " << sdig->GetId() << " " << sdig->GetAmp()/1.e+6 << endl;
188    }
189    }
190    }
191    
192    if(cbranch->GetEvent(0)) {
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){
197          cE->Fill(rp->GetEnergy());
198          cout << "RecPoint info " << rp->GetAbsId() << " " << rp->GetEnergy() << endl;
199    }
200    }
201    }
202    
203    */
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   
234 }