bug fix in the vertex selection
[u/mrichter/AliRoot.git] / PHOS / AliPHOSDApi0mip.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15 /* $Id$ */
16
17 // --
18 // --
19 // Implementation for TTree output in PHOS DA
20 // for calibrating energy by pi0 and MIP.
21 // --
22 // -- Author: Hisayuki Torii (Hiroshima Univ.)
23 // --
24
25 #include <stdio.h>
26 #include <iostream>
27 #include <math.h>
28 #include "AliLog.h"
29 #include "AliPHOSDApi0mip.h"
30 ClassImp(AliPHOSDApi0mip)
31 //----------------------------------------------------------------
32 AliPHOSDApi0mip::AliPHOSDApi0mip(int module,int iterid,char* fopt) :
33   TNamed(), fCreateTree(false), fCreateHist(false), fMod(0), fIterId(0),
34   fTFile(0), fTTree(0), fEvent(0), fEventClustered(0), fTime(0),
35   fH1Time(0), fH1DigitNum(0), fH1ClusterNum(0),
36   fH2EneDigitId(0), fH2MipDigitId(0), fH2Pi0DigitId(0), fH3Pi0AsymPt(0) {
37   // Constructor
38
39   char hname[1024], htitle[1024];
40   sprintf(hname,"AliPHOSDApi0mip_mod%d_iter%d",module,iterid);
41   SetName(hname);
42   sprintf(htitle,"PHOS MIP/pi0 Calibration DA for Module:%d Iteration:%d",module,iterid);
43   SetTitle(htitle);
44
45   fMod = module;
46   fIterId = iterid;
47   fEvent = 0;
48   fEventClustered = false;
49   fTime = 0;
50   fCreateTree = false;
51   fCreateHist = false;
52
53   char fname[1024];
54   sprintf(fname,"AliPHOSDApi0mip_mod%d.root",module);
55   fTFile = TFile::Open(fname,fopt);
56 }
57 //----------------------------------------------------------------
58 AliPHOSDApi0mip::AliPHOSDApi0mip(const AliPHOSDApi0mip& da):
59   TNamed(da), fCreateTree(false), fCreateHist(false), fMod(da.fMod), fIterId(da.fIterId),
60   fTFile(0), fTTree(0), fEvent(0), fEventClustered(0), fTime(0),
61   fH1Time(0), fH1DigitNum(0), fH1ClusterNum(0),
62   fH2EneDigitId(0), fH2MipDigitId(0), fH2Pi0DigitId(0), fH3Pi0AsymPt(0) {
63   // Copy Constructor
64
65   char fname[1024], hname[1024], htitle[1024];;
66   sprintf(fname,"%s.root",GetName());
67   fTFile = TFile::Open(fname,"RECREATE");
68   fEvent = 0;
69   fEventClustered = false;
70   fTime = 0;
71
72   if( da.fCreateHist ){
73     fCreateHist = true;
74     fH1Time = (TH1I*) da.fH1Time->Clone();
75     fH1DigitNum = (TH1F*) da.fH1DigitNum->Clone();
76     fH1ClusterNum = (TH1F*) da.fH1ClusterNum->Clone();
77     fH2EneDigitId = (TH2F*) da.fH2EneDigitId->Clone();
78     fH2MipDigitId = (TH2F*) da.fH2MipDigitId->Clone();
79     fH2Pi0DigitId = (TH2F*) da.fH2Pi0DigitId->Clone();
80     fH3Pi0AsymPt = (TH3F*) da.fH3Pi0AsymPt->Clone();
81   } else {
82     fCreateHist = false;
83     fH1Time = 0;
84     fH1DigitNum = 0;
85     fH1ClusterNum = 0;
86     fH2EneDigitId = 0;
87     fH2MipDigitId = 0;
88     fH2Pi0DigitId = 0;
89     fH3Pi0AsymPt = 0;
90   }
91   if( da.fCreateTree ){
92     // Create new ttree instead of copy.
93     fCreateTree = true;
94     sprintf(hname,"tevt_mod%d_iter%d",fMod,fIterId);
95     sprintf(htitle,"Calibration for Module:%d Iteration:%d",fMod,fIterId);
96     fTTree = new TTree(hname,htitle);
97     fTTree->Branch("AliPHOSDATreeEvent","AliPHOSDATreeEvent",&fEvent);
98   } else {
99     fCreateTree = false;
100     fTTree = 0;
101   }
102 }
103 //-------------------------------------------------------------------
104 AliPHOSDApi0mip& AliPHOSDApi0mip::operator= (const AliPHOSDApi0mip&)
105 {
106   // Operator= is not implemented yet
107   AliFatal("Operator= is not implemented");  
108   return *this;
109 }
110 //-------------------------------------------------------------------
111 AliPHOSDApi0mip::~AliPHOSDApi0mip(){
112   // Destructor
113
114   //std::cout<<" DEBUGDEBUG :: AliPHOSDApi0mip:: deconstructor "<<std::endl;
115   if( fCreateHist ){
116     fH1Time->Write(0,TObject::kOverwrite);
117     fH1DigitNum->Write(0,TObject::kOverwrite);
118     fH1ClusterNum->Write(0,TObject::kOverwrite);
119     fH2EneDigitId->Write(0,TObject::kOverwrite);
120     fH2MipDigitId->Write(0,TObject::kOverwrite);
121     fH2Pi0DigitId->Write(0,TObject::kOverwrite);
122     fH3Pi0AsymPt->Write(0,TObject::kOverwrite);
123     delete fH1Time;
124     delete fH1DigitNum;
125     delete fH1ClusterNum;
126     delete fH2EneDigitId;
127     delete fH2MipDigitId;
128     delete fH2Pi0DigitId;
129     delete fH3Pi0AsymPt;
130   }
131   if( fCreateTree ){
132     fTTree->Write();
133   }
134   if( fEvent ) delete fEvent;
135   fTFile->Close();
136   delete fTFile;
137 }
138 //-------------------------------------------------------------------
139 void AliPHOSDApi0mip::FillDigit(float adc,int row,int col){
140   // Put new digit information in event
141   //
142
143   if( fEvent==0 ) fEvent = new AliPHOSDATreeEvent;
144   float energy = adc * 0.005; // 5MeV/ch is fixed.
145   if( energy > 0.020 ){  // Threshold at 20MeV.
146     //std::cout<<" AliPHOSDApi0mip::AppendDigit() DEBUG:: (row,col)=("<<row<<","<<col<<") adc="<<adc<<" energy="<<energy<<std::endl;
147     fEvent->Fill(energy,row,col);
148   }
149 }
150 //-------------------------------------------------------------------
151 void AliPHOSDApi0mip::NewEvent(){
152   // Initiator must be called in the beginning of event
153   //
154
155   if(fEvent) fEvent->Reset();
156   fEventClustered = false;
157 }
158 //----------------------------------------------------------------
159 bool AliPHOSDApi0mip::CreateTree(){
160   if(! fTFile ) {
161     return false;
162   }
163   fTFile->cd();
164   char hname[1024], htitle[1024];
165   sprintf(hname,"trevt_mod%d_iter%d",fMod,fIterId);
166   sprintf(htitle,"Calibration for Module:%d Iteration:%d",fMod,fIterId);
167   fTTree = new TTree(hname,htitle);
168   fTTree->Branch("AliPHOSDATreeEvent","AliPHOSDATreeEvent",&fEvent);
169   fCreateTree = true;
170   return true;
171 }
172 //----------------------------------------------------------------
173 bool AliPHOSDApi0mip::CreateHist(){
174   // Create histogram routine called by the constructor only.
175
176   char hname[1024], htitle[1024];
177
178   sprintf(hname,"h1_time_mod%d_iter%d",fMod,fIterId);
179   sprintf(htitle,"Time : Mod:%d Iter:%d",fMod,fIterId);
180   fH1Time = (TH1I*) gDirectory->Get(hname);
181   if( fH1Time>0 ){
182     std::cout<<" AliPHOSDApi0mip:Warning!! Output object already exist : "<<fH1Time->GetName()<<std::endl;
183     return false;
184   }
185   fH1Time = new TH1I(hname,htitle,2,0,2);
186
187   sprintf(hname,"h1_digitnum_mod%d_iter%d",fMod,fIterId);
188   sprintf(htitle,"Number of Digits : Mod:%d Iter:%d",fMod,fIterId);
189   fH1DigitNum = new TH1F(hname,htitle,100,0,100);
190
191   sprintf(hname,"h1_clusternum_mod%d_iter%d",fMod,fIterId);
192   sprintf(htitle,"Number of Clusters : Mod:%d Iter:%d",fMod,fIterId);
193   fH1ClusterNum = new TH1F(hname,htitle,100,0,100);
194
195   sprintf(hname,"h2_pi0digitid_mod%d_iter%d",fMod,fIterId);
196   sprintf(htitle,"PHOS pi0 mass vs Digit Id : Mod:%d Iter:%d",fMod,fIterId);
197   fH2Pi0DigitId = new TH2F(hname,htitle,17920,0,17920,120,0,0.3);
198       
199   sprintf(hname,"h2_mipdigitid_mod%d_iter%d",fMod,fIterId);
200   sprintf(htitle,"PHOS MIP vs Digit Id : Mod:%d Iter:%d",fMod,fIterId);
201   fH2MipDigitId = new TH2F(hname,htitle,17920,0,17920,100,0-0.0025,0.5-0.0025);
202   //fH2MipDigitId = new TH2F(hname,htitle,17920,0,17920,50,0.5,1.5);
203
204   sprintf(hname,"h2_enedigitid_mod%d_iter%d",fMod,fIterId);
205   sprintf(htitle,"PHOS MIP vs Digit Id : Mod:%d Iter:%d",fMod,fIterId);
206   //fH2EneDigitId = new TH2F(hname,htitle,17920,0,17920,100,0-0.0025,0.5-0.0025);
207   fH2EneDigitId = new TH2F(hname,htitle,17920,0,17920,50,0,5.0);
208
209   sprintf(hname,"h3_pi0asympt_mod%d_iter%d",fMod,fIterId);
210   sprintf(htitle,"PHOS pi0 mass vs pT vs Asym : Mod:%d Iter:%d",fMod,fIterId);
211   fH3Pi0AsymPt = new TH3F(hname,htitle,20,0,1,20,0,10,200,0,1);
212
213   fCreateHist = true;
214   return true;
215 }
216 //-------------------------------------------------------------------
217 void AliPHOSDApi0mip::FillTree(AliPHOSDATreeEvent* event){
218   // Fill all information (AliPHOSDATreeEvent) into TTree
219   // if event==NULL, the obtained information (=fEvent) is used.
220
221   if(! fCreateTree ) CreateTree();
222   if( event == 0 ){
223     if( fEvent == 0 ){
224       fEvent = new AliPHOSDATreeEvent();
225     }
226     // No need to make clustering for tree saving.
227     //if(! fEventClustered ){
228     //fEvent->ExecuteClustering();
229     //fEventClustered = true;
230     //}
231     fEvent->SetTime(fTime);
232     fTTree->Fill();
233   } else {
234     if( fEvent ) delete fEvent;
235     fEvent = event;
236     fTTree->Fill();
237     fEvent = 0;
238   }
239 }
240 //-------------------------------------------------------------------
241 void AliPHOSDApi0mip::FillHist(AliPHOSDATreeEvent* event){
242   // Fill all information (AliPHOSDATreeEvent) into histograms
243   // if event==NULL, the obtained information (=fEvent) is used.
244
245   if(! fCreateHist ) CreateHist();
246   if( event == 0 ){
247     if( fEvent==0 ) fEvent = new AliPHOSDATreeEvent;
248     if(! fEventClustered ){
249       fEvent->ExecuteClustering();
250       fEventClustered = true;
251     }
252     fEvent->SetTime(fTime);
253     event = fEvent;
254   }
255
256   // -- Constant
257   int nclt1, nclt2;
258   int ndigits;
259   int ndigitsall = 0;
260   AliPHOSDATreeDigit digit1, digit2;
261   const float kDistanceFromVertex = 460.;
262   float dist1, dist2, px, py;
263   float mass, cosine, pt, asym;
264
265   // -- Filling into Time Histogram
266   double start = fH1Time->GetBinContent(1);
267   double stop = fH1Time->GetBinContent(2);
268   double currenttime = (double) (event->GetTime());
269   if( start == 0 && stop == 0 ) start = stop = currenttime;
270   if( currenttime < start ) start = currenttime;
271   if( currenttime > stop ) stop = currenttime;
272   fH1Time->SetBinContent(1,start);
273   fH1Time->SetBinContent(2,stop);
274
275   // -- Start looping clusters
276   nclt1 = event->GetNClusters();
277   fH1ClusterNum->Fill(nclt1);
278   //std::cout<<" AliPHOSDApi0mip::FillHisto() DEBUG:: Number of clusters: "<<nclt1<<std::endl;
279   while( nclt1-- ){
280     AliPHOSDATreeCluster& clt1 = event->GetCluster(nclt1);
281     ndigits = clt1.GetNDigits();
282     while( ndigits-- ){
283       AliPHOSDATreeDigit& digit = clt1.GetDigit(ndigits);
284       fH2EneDigitId->Fill(digit.GetDigitId(),digit.GetEnergy());
285       ndigitsall++;
286     }
287     digit1 = clt1.GetMaxDigit();
288     if( clt1.GetEnergy() > 0 && clt1.GetNDigits() >= 2 ){
289       fH2MipDigitId->Fill(digit1.GetDigitId(),clt1.GetEnergy());
290     }
291     if( clt1.GetEnergy() > 0.2 ){
292       nclt2 = nclt1;
293       while( nclt2-- ){
294         AliPHOSDATreeCluster& clt2 = event->GetCluster(nclt2);
295         digit2 = clt2.GetMaxDigit();
296         if( clt2.GetEnergy() > 0.2 ){
297           // 
298           // Following approximation is applied in order to avoid many calling "sqrt" function.
299           // sqrt( 1 + x^2 ) =~ 1 + 0.5x^2 (for x<<1)
300           //
301           //cosine = ((clt1.GetRow()-31.5)*(clt2.GetRow()-31.5)*4.84 + (clt1.GetCol()-27.5)*(clt2.GetCol()-27.5)*4.84 + kDistanceFromVertex*kDistanceFromVertex)
302           //        / ( kDistanceFromVertex + (clt1.GetRow()-31.5)*(clt1.GetRow()-31.5)*2.42/kDistanceFromVertex + (clt1.GetCol()-27.5)*(clt1.GetCol()-27.5)*2.42/kDistanceFromVertex )
303           //        / ( kDistanceFromVertex + (clt2.GetRow()-31.5)*(clt2.GetRow()-31.5)*2.42/kDistanceFromVertex + (clt2.GetCol()-27.5)*(clt2.GetCol()-27.5)*2.42/kDistanceFromVertex );
304           //
305           
306
307           dist1 = sqrt( kDistanceFromVertex*kDistanceFromVertex + (clt1.GetRow()-31.5)*(clt1.GetRow()-31.5)*4.84 + (clt1.GetCol()-27.5)*(clt1.GetCol()-27.5)*4.84 );
308           dist2 = sqrt( kDistanceFromVertex*kDistanceFromVertex + (clt2.GetRow()-31.5)*(clt2.GetRow()-31.5)*4.84 + (clt2.GetCol()-27.5)*(clt2.GetCol()-27.5)*4.84 );
309           cosine = ((clt1.GetRow()-31.5)*(clt2.GetRow()-31.5)*4.84 + (clt1.GetCol()-27.5)*(clt2.GetCol()-27.5)*4.84 + kDistanceFromVertex*kDistanceFromVertex)
310             / dist1 / dist2;
311           mass = sqrt( 2.0 * clt1.GetEnergy() * clt2.GetEnergy() * (1 - cosine) );
312           //std::cout<<" DEBUG: dist1:"<<dist1<<" dist2:"<<dist2<<" cosine:"<<cosine<<" c1:"<<clt1.GetEnergy()<<" c2:"<<clt2.GetEnergy()<<" mass:"<<mass<<std::endl;
313           //std::cout<<"      : clt1(row,col)=("<<clt1.GetRow()<<","<<clt1.GetCol()<<") clt2(row,col)=("<<clt2.GetRow()<<","<<clt2.GetCol()<<") "<<std::endl;
314           asym = fabs(clt1.GetEnergy() - clt2.GetEnergy())/(clt1.GetEnergy() + clt2.GetEnergy());
315           py = clt1.GetEnergy() * kDistanceFromVertex / dist1 + clt2.GetEnergy() * kDistanceFromVertex / dist2;
316           px = clt1.GetEnergy() * (clt1.GetCol()-27.5)*2.42/dist1 + clt2.GetEnergy() * (clt2.GetCol()-27.5)*2.42/ dist2;
317           pt = sqrt(px*px+py*py);
318           fH3Pi0AsymPt->Fill(asym,pt,mass);
319           if( digit1.IsValid() && digit1.IsValid() &&
320               digit1.GetDigitId() != digit2.GetDigitId() &&
321               //pt > 1.5 && clt1.GetEnergy() > 0.3 && clt2.GetEnergy() > 0.3 ){
322               //pt > 1.3 && clt1.GetEnergy() > 0.25 && clt2.GetEnergy() > 0.25 ){
323               pt > 1.0 && clt1.GetEnergy() > 0.2 && clt2.GetEnergy() > 0.2 ){
324             fH2Pi0DigitId->Fill(digit1.GetDigitId(),mass);
325             fH2Pi0DigitId->Fill(digit2.GetDigitId(),mass);
326           }
327         } // End of energy cut
328       } // End of second cluster loop (clt2)
329     } // End of energy cut
330   } // End of first cluster loop (clt1)
331   fH1DigitNum->Fill(ndigitsall);
332   //
333 }
334 //-------------------------------------------------------------------
335 void AliPHOSDApi0mip::Print(Option_t *option) const
336 {
337   // Print Out
338
339   //fTFile->ls();
340   //fTTree->Print();
341   if( fEvent ){
342     //fEvent->ExecuteClustering();
343     fEvent->Print(option);
344   }
345 }
346 //-------------------------------------------------------------------