MIP calculation added
[u/mrichter/AliRoot.git] / VZERO / AliVZEROReconstructor.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 /* $Id$ */
17
18 ///////////////////////////////////////////////////////////////////////////////
19 ///                                                                          //
20 /// class for VZERO reconstruction                                           //
21 ///                                                                          //
22 ///////////////////////////////////////////////////////////////////////////////
23
24 #include "AliRunLoader.h"
25 #include "AliRawReader.h"
26 #include "AliVZEROReconstructor.h"
27 #include "AliVZERORawStream.h"
28 #include "AliESDEvent.h"
29 #include "AliVZEROTriggerMask.h"
30
31 ClassImp(AliVZEROReconstructor)
32
33 //_____________________________________________________________________________
34 AliVZEROReconstructor:: AliVZEROReconstructor(): AliReconstructor(),
35    fESDVZERO(0x0),
36    fESD(0x0),
37    fCalibData(GetCalibData())
38 {
39   // Default constructor  
40   // Get calibration data
41   
42   // fCalibData = GetCalibData(); 
43
44 }
45
46
47 //_____________________________________________________________________________
48 AliVZEROReconstructor& AliVZEROReconstructor::operator = 
49   (const AliVZEROReconstructor& /*reconstructor*/)
50 {
51 // assignment operator
52
53   Fatal("operator =", "assignment operator not implemented");
54   return *this;
55 }
56
57 //_____________________________________________________________________________
58 AliVZEROReconstructor::~AliVZEROReconstructor()
59 {
60 // destructor
61
62    delete fESDVZERO;    
63 }
64
65 //_____________________________________________________________________________
66 void AliVZEROReconstructor::Init()
67 {
68 // initializer
69
70   fESDVZERO  = new AliESDVZERO;
71 }
72
73 //______________________________________________________________________
74 void AliVZEROReconstructor::ConvertDigits(AliRawReader* rawReader, TTree* digitsTree) const
75 {
76 // converts RAW to digits - pedestal is subtracted 
77
78   if (!digitsTree) {
79     AliError("No digits tree!");
80     return;
81   }
82
83   TClonesArray* digitsArray = new TClonesArray("AliVZEROdigit");
84   digitsTree->Branch("VZERODigit", &digitsArray);
85
86   rawReader->Reset();
87   AliVZERORawStream rawStream(rawReader);
88   if (rawStream.Next()) {  
89      Int_t ADC_max[64], adc[64], time[64];   
90      for(Int_t i=0; i<64; i++) {
91          // Search for the maximum charge in the train of 21 LHC clocks 
92          // regardless of the integrator which has been operated:
93          ADC_max[i] = 0;
94          for(Int_t iClock=0; iClock<21; iClock++){
95              if((Int_t)rawStream.GetPedestal(i,iClock) > ADC_max[i])  
96                 {ADC_max[i]=(Int_t)rawStream.GetPedestal(i,iClock);}
97          }
98          // Convert i (FEE channel numbering) to j (aliroot channel numbering)
99          Int_t j =  rawStream.GetOfflineChannel(i);
100          adc[j]  =  ADC_max[i];
101          time[j] =  rawStream.GetTime(i); 
102      }  
103      // Channels(aliroot numbering) will be ordered in the tree
104      for(Int_t iChannel = 0; iChannel < 64; iChannel++) {
105          new ((*digitsArray)[digitsArray->GetEntriesFast()])
106              AliVZEROdigit(iChannel,adc[iChannel],time[iChannel]);
107      }
108   }
109
110   digitsTree->Fill();
111 }
112
113 //______________________________________________________________________
114 void AliVZEROReconstructor::FillESD(TTree* digitsTree, TTree* /*clustersTree*/,
115                                     AliESDEvent* esd) const
116 {
117 // fills multiplicities to the ESD
118
119   if (!digitsTree) {
120     AliError("No digits tree!");
121     return;
122   }
123
124   TClonesArray* digitsArray = NULL;
125   TBranch* digitBranch = digitsTree->GetBranch("VZERODigit");
126   digitBranch->SetAddress(&digitsArray);
127
128   Short_t Multiplicity[64];
129   Float_t   mult[64];  
130   Short_t    adc[64]; 
131   Short_t   time[64]; 
132   for (Int_t i=0; i<64; i++){
133        adc[i] = 0;
134        mult[i]= 0.0;
135   }
136      
137   // loop over VZERO entries to get multiplicity
138   Int_t nEntries = (Int_t)digitsTree->GetEntries();
139   for (Int_t e=0; e<nEntries; e++) {
140     digitsTree->GetEvent(e);
141
142     Int_t nDigits = digitsArray->GetEntriesFast();
143     
144     for (Int_t d=0; d<nDigits; d++) {    
145         AliVZEROdigit* digit = (AliVZEROdigit*)digitsArray->At(d);      
146         Int_t  pmNumber      = digit->PMNumber(); 
147         // Pedestal retrieval and suppression: 
148         Int_t  pedestal      = int(fCalibData->GetPedestal(d));
149         adc[pmNumber]  = (Short_t) digit->ADC() - pedestal; 
150         time[pmNumber] = (Short_t) digit->Time();
151         // printf("PM = %d,  MIP per ADC channel = %f \n",pmNumber, GetMIP(pmNumber));
152         // cut of ADC at 1MIP/2 
153         if (adc[pmNumber] > (int(1.0/GetMIP(pmNumber)) /2) ) 
154             mult[pmNumber] += float(adc[pmNumber])*GetMIP(pmNumber);
155     } // end of loop over digits
156   } // end of loop over events in digits tree
157   
158   for (Int_t j=0; j<64; j++) Multiplicity[j] = short(mult[j]+0.5);       
159   fESDVZERO->SetMultiplicity(Multiplicity);
160   fESDVZERO->SetADC(adc);
161   fESDVZERO->SetTime(time);
162
163   // now get the trigger mask
164
165   AliVZEROTriggerMask *TriggerMask = new AliVZEROTriggerMask();
166   TriggerMask->SetAdcThreshold(20.0/2.0);
167   TriggerMask->SetTimeWindowWidthBBA(50);
168   TriggerMask->SetTimeWindowWidthBGA(20);
169   TriggerMask->SetTimeWindowWidthBBC(50);
170   TriggerMask->SetTimeWindowWidthBGC(20);
171   TriggerMask->FillMasks(digitsTree,digitsArray);
172
173   fESDVZERO->SetBBtriggerV0A(TriggerMask->GetBBtriggerV0A());
174   fESDVZERO->SetBGtriggerV0A(TriggerMask->GetBGtriggerV0A());
175   fESDVZERO->SetBBtriggerV0C(TriggerMask->GetBBtriggerV0C());
176   fESDVZERO->SetBGtriggerV0C(TriggerMask->GetBGtriggerV0C());
177   
178   if (esd) { 
179      AliDebug(1, Form("Writing VZERO data to ESD tree"));
180      esd->SetVZEROData(fESDVZERO);
181   }
182 }
183
184 //_____________________________________________________________________________
185 AliCDBStorage* AliVZEROReconstructor::SetStorage(const char *uri) 
186 {
187 // Sets the storage  
188
189   Bool_t deleteManager = kFALSE;
190   
191   AliCDBManager *manager = AliCDBManager::Instance();
192   AliCDBStorage *defstorage = manager->GetDefaultStorage();
193   
194   if(!defstorage || !(defstorage->Contains("VZERO"))){ 
195      AliWarning("No default storage set or default storage doesn't contain VZERO!");
196      manager->SetDefaultStorage(uri);
197      deleteManager = kTRUE;
198   }
199  
200   AliCDBStorage *storage = manager->GetDefaultStorage();
201
202   if(deleteManager){
203      AliCDBManager::Instance()->UnsetDefaultStorage();
204      defstorage = 0;   // the storage is killed by AliCDBManager::Instance()->Destroy()
205   }
206
207   return storage; 
208 }
209
210 //_____________________________________________________________________________
211 AliVZEROCalibData* AliVZEROReconstructor::GetCalibData() const
212 {
213   // Gets calibration object for VZERO set
214
215   AliCDBManager *man = AliCDBManager::Instance();
216
217   AliCDBEntry *entry=0;
218
219   entry = man->Get("VZERO/Calib/Data");
220
221 //   if(!entry){
222 //     AliWarning("Load of calibration data from default storage failed!");
223 //     AliWarning("Calibration data will be loaded from local storage ($ALICE_ROOT)");
224 //     Int_t runNumber = man->GetRun();
225 //     entry = man->GetStorage("local://$ALICE_ROOT")
226 //       ->Get("VZERO/Calib/Data",runNumber);
227 //      
228 //   }
229
230   // Retrieval of data in directory VZERO/Calib/Data:
231
232   AliVZEROCalibData *calibdata = 0;
233
234   if (entry) calibdata = (AliVZEROCalibData*) entry->GetObject();
235   if (!calibdata)  AliFatal("No calibration data from calibration database !");
236
237   return calibdata;
238 }
239
240 //_____________________________________________________________________________
241 Float_t AliVZEROReconstructor::GetMIP(Int_t channel) const {
242
243 // Computes the MIP conversion factor - MIP per ADC channel - 
244 // Argument passed is the PM number (aliroot numbering)
245
246   Float_t P0[64] = {
247   7.094891, 7.124938, 7.089708, 7.098169, 7.094482, 7.147250, 7.170978, 7.183392, 
248   7.145760, 7.148096, 7.153840, 7.143544, 7.186069, 7.194580, 7.203516, 7.195176, 
249   7.188333, 7.198607, 7.209412, 7.226565, 7.221695, 7.205132, 7.191238, 7.227724, 
250   7.232810, 7.252655, 7.230309, 7.273518, 7.273518, 7.242969, 7.252859, 7.252655, 
251   7.026802, 7.079913, 7.134147, 7.092387, 7.079561, 7.072848, 7.123192, 7.003141, 
252   7.024667, 7.124784, 7.123442, 7.129744, 7.110671, 7.143031, 7.139439, 7.178109, 
253   7.247803, 7.139396, 7.293809, 7.094454, 6.992198, 7.206448, 7.244765, 7.056197, 
254   7.263595, 7.138569, 7.089582, 7.215683, 7.266183, 7.165123, 7.243276, 7.235135 };
255   Float_t P1[64] = {
256   0.135569, 0.146405, 0.142425, 0.144278, 0.142307, 0.141648, 0.128477, 0.138239, 
257   0.144173, 0.143419, 0.143572, 0.144482, 0.138024, 0.136542, 0.135955, 0.138537, 
258   0.148521, 0.141999, 0.139627, 0.130014, 0.134970, 0.135635, 0.139094, 0.140634, 
259   0.137971, 0.142080, 0.142793, 0.142778, 0.142778, 0.146045, 0.139133, 0.142080, 
260   0.144121, 0.142311, 0.136564, 0.142686, 0.138792, 0.166285, 0.136387, 0.155391, 
261   0.176082, 0.140408, 0.164738, 0.144270, 0.142766, 0.147486, 0.141951, 0.138012, 
262   0.132394, 0.142849, 0.140477, 0.144592, 0.141558, 0.157646, 0.143758, 0.173385, 
263   0.146489, 0.143279, 0.145230, 0.147203, 0.147333, 0.144979, 0.148597, 0.138985 };
264
265 // High Voltage retrieval from Calibration Data Base:  
266   Float_t  HV = fCalibData->GetMeanHV(channel);  
267   Float_t MIP = 0.5/TMath::Exp((TMath::Log(HV) - P0[channel] )/P1[channel]);
268   return MIP; 
269
270 }