LHAPDF 5.5.1 and 5.9.1
[u/mrichter/AliRoot.git] / FIT / AliFITReconstructor.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
16 /*********************************************************************
17  *  FIT reconstruction and filling ESD
18  *  - reconstruct mean time (interation time) 
19  *  - vertex position
20  *  -  multiplicity
21  ********************************************************************/
22
23 #include <AliESDEvent.h>
24 #include "AliLog.h"
25 #include "AliFITRecPoint.h"
26 #include "AliFITDigit.h"
27 #include "AliFITReconstructor.h"
28 #include "AliESDFIT.h"
29 #include "AliLog.h"
30 #include "AliESDRun.h"
31 #include "AliFITRawReader.h"
32
33 #include <TArrayI.h>
34 #include <TGraph.h>
35 #include <TMath.h>
36 #include <Riostream.h>
37
38 using std::cout;
39 using std::endl;
40
41 ClassImp(AliFITReconstructor)
42
43 AliFITReconstructor:: AliFITReconstructor(): AliReconstructor(),
44                                              fESD(NULL),
45                                              fESDFIT(NULL),
46                                              fDigits(NULL)
47
48 {
49  
50   //constructor
51  
52   
53    fESDFIT  = new AliESDFIT();
54
55  
56 }
57 //_____________________________________________________________________________
58 void AliFITReconstructor::Init()
59 {
60 // initializer
61
62   fESDFIT  = new AliESDFIT;
63  }
64
65 //______________________________________________________________________
66 void AliFITReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digitsTree) const
67 {
68 // converts RAW to digits 
69   cout<<" AliFITReconstructor::ConvertDigits "<<rawReader<<" "<<digitsTree<<endl;
70   if (!digitsTree) {
71     AliError("No digits tree!");
72     return;
73   }
74
75   fDigits = new TClonesArray ("AliFITDigit", 100);
76   digitsTree->Branch("FIT", &fDigits);
77
78   AliFITRawReader myrawreader(rawReader);
79   if (myrawreader.Next()) 
80     {
81       Int_t allData[500];
82       for (Int_t i=0; i<500; i++)  allData[i]=0;
83       for (Int_t i=0; i<500; i++) 
84         if(myrawreader.GetData(i)>0)  { 
85           allData[i]=myrawreader.GetData(i);
86         }
87    
88       Int_t timeCFD, timeLED, timeQT1, timeQT0;
89       for (Int_t ipmt=0; ipmt<160; ipmt++) {
90         if(allData[ipmt]>0) {
91           timeCFD = allData[ipmt];
92           timeLED = allData[ipmt];
93           timeQT0= allData[ipmt+160];
94           timeQT1 = allData[ipmt+320];
95           //add digit
96           new((*fDigits)[fDigits->GetEntriesFast()] )
97             AliFITDigit( ipmt,timeCFD, timeLED, timeQT0, timeQT1, 0);
98         }
99       }
100       
101       
102       
103       digitsTree->Fill(); 
104     }
105
106   
107 }
108
109  //____________________________________________________________
110 void AliFITReconstructor::FillESD(TTree *digitsTree, TTree * /*clustersTree*/, AliESDEvent *pESD) const
111 {
112   
113   /***************************************************
114   Resonstruct digits to vertex position
115   ****************************************************/
116   
117   AliDebug(1,Form("Start FillESD FIT"));
118   if(!pESD) {
119     AliError("No ESD Event");
120     return;
121   }
122   
123   Float_t channelWidth = 24.4;  
124   Float_t c = 0.0299792458; // cm/ps
125   Float_t currentVertex, shift=0;
126   Int_t ncont=-1;
127    const AliESDVertex* vertex = pESD->GetPrimaryVertex();
128   if (!vertex)        vertex = pESD->GetPrimaryVertexSPD();
129   if (!vertex)        vertex = pESD->GetPrimaryVertexTPC();
130   if (!vertex)        vertex = pESD->GetVertex();
131   
132   if (vertex) {
133     AliDebug(2, Form("Got %s (%s) from ESD: %f", 
134                      vertex->GetName(), vertex->GetTitle(), vertex->GetZ()));
135     currentVertex = vertex->GetZ();
136     
137     ncont = vertex->GetNContributors();
138     if(ncont>0 ) {
139       shift = currentVertex/c;
140     }
141   } //vertex
142
143  
144   // FIT digits reconstruction
145   
146   if (!digitsTree) {
147     AliError("No digits tree!");
148     return;
149   }
150   
151   fDigits = new TClonesArray("AliFITDigit", 100);
152   // digitsTree->Print();
153   TBranch *brDigits=digitsTree->GetBranch("FIT");
154   if (brDigits) 
155     brDigits->SetAddress(&fDigits);
156   else {
157     AliError("no FIT branch");
158     return;
159   }
160   const Float_t max_time = 1e7;
161   Int_t pmt=-1;
162   Float_t  time[160],amp[160];  
163   for(Int_t i=0; i<160; i++)   {  time[i]=max_time; amp[i]=0;}
164   
165   Int_t nEntries = (Int_t)digitsTree->GetEntries();
166   for (Int_t iev=0; iev<nEntries; iev++) {
167     digitsTree->GetEvent(iev,1);
168     brDigits->GetEvent(iev);
169     Int_t nDigits = fDigits->GetEntriesFast(); 
170     for (Int_t dig=0; dig<nDigits; dig++) {    
171       AliFITDigit* digit = (AliFITDigit*) fDigits->At(dig);      
172       pmt = digit->NPMT();
173       time[pmt] = Float_t (digit->TimeCFD() );
174       amp[pmt] = Float_t (digit->TimeQT1() - digit->TimeQT0() );
175     } 
176     fESDFIT->SetFITtime(time);         // best TOF on each PMT 
177     fESDFIT->SetFITamplitude(amp);     // number of particles(MIPs) on each 
178     
179     Float_t firsttime[3] = {max_time,max_time,max_time};
180     
181     Float_t vertexFIT = 9999999;
182     for (Int_t ipmt=0; ipmt<80; ipmt++)//timeC
183           if(time[ipmt]<firsttime[2]) firsttime[2]=time[ipmt]; 
184     
185     for ( Int_t ipmt=80; ipmt<160; ipmt++) 
186           if(time[ipmt]<firsttime[1]) firsttime[1]=time[ipmt]; 
187     
188     
189     AliDebug(5,Form("1stimeA %f , 1sttimeC %f  ",
190                     firsttime[1],
191                     firsttime[2] ) );
192     if (firsttime[1]<max_time && firsttime[2]<max_time)  {
193         firsttime[0] =  channelWidth *(firsttime[1] + firsttime[2])/2;
194         vertexFIT =  c*channelWidth*(firsttime[1] - firsttime[2])/2;
195     }
196       if (firsttime[1]<max_time) 
197         firsttime[1] = firsttime[1] * channelWidth + shift; 
198       
199       if (firsttime[2]<max_time) 
200         firsttime[2] = firsttime[2] * channelWidth - shift; 
201       
202       
203       
204       for(Int_t i=0; i<3; i++) {
205         fESDFIT->SetFITT0(i,firsttime[i]);   // interaction time (ns) 
206       }
207       fESDFIT->SetFITzVertex(vertexFIT); //vertex Z position 
208       
209       
210       AliDebug(1,Form("FIT: SPDshift %f Vertex %f  FITsignal %f ps FITA %f ps FITC %f ps \n",shift, vertexFIT, firsttime[0], firsttime[1],firsttime[2]));
211   }
212   
213    }