Update of calibration task and replace AliTRDrunCalib.C by AddTaskTRDCalib.C (Raphaelle)
[u/mrichter/AliRoot.git] / TRD / Macros / AliTRDCreate.C
1 #if !defined( __CINT__) || defined(__MAKECINT__)
2
3 #include <iostream>
4 #include <TRandom.h>
5 #include <TSystem.h>
6 #include <TDatime.h>
7
8 #include <AliCDBManager.h>
9 #include <AliCDBStorage.h>
10 #include <AliCDBEntry.h>
11 #include <AliCDBMetaData.h>
12 #include <AliGeometry.h>
13 #include <AliPID.h>
14
15 #include "../TRD/AliTRDgeometry.h"
16
17 #include "../TRD/Cal/AliTRDCalROC.h"
18 #include "../TRD/Cal/AliTRDCalPad.h"
19 #include "../TRD/Cal/AliTRDCalDet.h"
20
21 #endif
22
23 // run number for the dummy file
24 const Int_t gkDummyRun = 0;
25 AliCDBStorage* gStorLoc = 0;
26 Double_t t0minimum[540];
27 Double_t gainav[540];
28 Double_t vav[540];
29
30 TObject* CreatePadObjectT0Random(const char* shortName, const char* description, Float_t mean, Float_t sigma);
31 TObject* CreateDetT0Object(const char* shortName, const char* description);
32 Double_t funcpoly16(Double_t* x);
33 Double_t funcpoly12(Double_t* x);
34 Double_t funcpoly144(Double_t* x);
35 Double_t funcpoly16v(Double_t* x);
36 Double_t funcpoly12v(Double_t* x);
37 Double_t funcpoly144v(Double_t* x);
38 TObject* CreatePadObject(const char* shortName, const char* description, Float_t value);
39 TObject* CreatePadObjectRandom(const char* shortName, const char* description, Float_t mean, Float_t sigma);
40 TObject* CreatePadObjectRandomv(const char* shortName, const char* description, Float_t mean, Float_t sigma);
41 TObject* CreatePadObjectbridge(const char* shortName, const char* description);
42 TObject* CreatePadObjectbridgev(const char* shortName, const char* description);
43 TObject* CreateDetObject(const char* shortName, const char* description, Float_t value);
44 TObject* CreateDetObjectRandom(const char* shortName, const char* description, Float_t mean, Float_t sigma);
45 TObject* CreateDetObjectRandomg(const char* shortName, const char* description, Float_t mean, Float_t sigma);
46 TObject* CreateDetObjectRandomv(const char* shortName, const char* description, Float_t mean, Float_t sigma);
47 AliCDBMetaData* CreateMetaObject(const char* objectClassName);
48 void StoreObject(const char* cdbPath, TObject* object, AliCDBMetaData* metaData);
49 void AliTRDCreate(Bool_t residual = kFALSE);
50
51
52 //________________________________________________________________________________________________
53 TObject* CreatePadObjectT0Random(const char* shortName, const char* description, Float_t mean, Float_t sigma)
54 {
55   Double_t min = 0.0;
56   AliTRDCalPad *calPad = new AliTRDCalPad(shortName, description);
57   for (Int_t det=0; det<AliTRDgeometry::kNdet; ++det)
58     {
59       AliTRDCalROC *calROC = calPad->GetCalROC(det);
60       Double_t value[2400];
61       for (Int_t channel=0; channel<calROC->GetNchannels(); ++channel){
62         value[channel] = gRandom->Gaus(mean,sigma);
63         if(channel == 0) {
64           min = value[0];
65           //printf("value[0] %f and min %f et gRandom %f\n",value[0],min, gRandom->Gaus(mean,sigma));
66         }
67         if(min > value[channel]) min = value[channel];
68         //calROC->SetValue(channel, TMath::Abs(value));
69       }
70       t0minimum[det] = min;
71       for (Int_t channel=0; channel<calROC->GetNchannels(); ++channel){
72         //printf("min value for the det %d is %f and the channel is %f, we put %f\n",det,t0minimum[det],value[channel],(value[channel]-t0minimum[det]));
73         calROC->SetValue(channel, (value[channel]-t0minimum[det]));
74       }
75     }
76   
77   return calPad;
78 }
79 //___________________________________________________________________________________________________
80 TObject* CreateDetT0Object(const char* shortName, const char* description)
81 {
82   AliTRDCalDet *object = new AliTRDCalDet(shortName, description);
83   for (Int_t det=0; det<AliTRDgeometry::kNdet; ++det)
84     object->SetValue(det, t0minimum[det]);
85   return object;
86 }
87 //__________________________________________________________________________________________________
88 Double_t funcpoly16(Double_t* x){
89
90   //0.1
91
92   Double_t par[3] = {-0.002678571429,0.040178571,1};
93
94   // sum landau + gaus with identical mean
95
96   Double_t valLandau = par[0]*x[0]*x[0]+par[1]*x[0]+par[2];
97
98   return valLandau;
99 }
100 //_____________________________________________________________________________________________________
101 Double_t funcpoly12(Double_t* x){
102
103   // 0.1
104
105   Double_t par[3] = {-0.005,0.055,1};
106
107   // sum landau + gaus with identical mean
108
109   Double_t valLandau = par[0]*x[0]*x[0]+par[1]*x[0]+par[2];
110
111   return valLandau;
112 }
113 //______________________________________________________________________________________________________
114 Double_t funcpoly144(Double_t* x){
115
116   //0.1
117
118   Double_t par[3] = {-0.00001956181536,0.00279339596,1};
119
120   // sum landau + gaus with identical mean
121
122   Double_t valLandau = par[0]*x[0]*x[0]+par[1]*x[0]+par[2];
123
124   return valLandau;
125 }
126
127 //__________________________________________________________________________________________________
128 Double_t funcpoly16v(Double_t* x){
129
130   //0.05 
131
132   Double_t par[3] = {-0.0008928571429,0.013392857,1};
133
134   // sum landau + gaus with identical mean
135
136   Double_t valLandau = par[0]*x[0]*x[0]+par[1]*x[0]+par[2];
137
138   return valLandau;
139 }
140 //_____________________________________________________________________________________________________
141 Double_t funcpoly12v(Double_t* x){
142
143   //0.05
144
145   Double_t par[3] = {-0.001666667,0.018333,1};
146
147   // sum landau + gaus with identical mean
148
149   Double_t valLandau = par[0]*x[0]*x[0]+par[1]*x[0]+par[2];
150
151   return valLandau;
152 }
153 //___________________________________________________________________________________________________
154 Double_t funcpoly144v(Double_t* x){
155
156   //0.05
157
158   Double_t par[3] = {-0.000009780907668,0.001398669797,1};
159
160   // sum landau + gaus with identical mean
161
162   Double_t valLandau = par[0]*x[0]*x[0]+par[1]*x[0]+par[2];
163
164   return valLandau;
165 }
166 //__________________________________________________________________________________________________
167 TObject* CreatePadObject(const char* shortName, const char* description, Float_t value)
168 {
169   AliTRDCalPad *calPad = new AliTRDCalPad(shortName, description);
170   for (Int_t det=0; det<AliTRDgeometry::kNdet; ++det)
171     {
172       AliTRDCalROC *calROC = calPad->GetCalROC(det);
173       for (Int_t channel=0; channel<calROC->GetNchannels(); ++channel){
174         calROC->SetValue(channel, value);
175       }
176     }
177   return calPad;
178 }
179 //___________________________________________________________________________________________________
180 TObject* CreatePadObjectRandom(const char* shortName, const char* description, Float_t mean, Float_t sigma)
181 {
182   AliTRDCalPad *calPad = new AliTRDCalPad(shortName, description);
183   for (Int_t det=0; det<AliTRDgeometry::kNdet; ++det)
184     {
185     
186       AliTRDCalROC *calROC = calPad->GetCalROC(det);
187       for (Int_t channel=0; channel<calROC->GetNchannels(); ++channel){
188         Double_t value = gRandom->Gaus(mean,sigma);
189         //cout << "value: " << value << endl;
190         //if(value < 0) calROC->SetValue(channel, 0);
191         calROC->SetValue(channel, value);
192       }
193     }
194   
195   return calPad;
196 }
197 //_______________________________________________________________________________________
198 TObject* CreatePadObjectRandomv(const char* shortName, const char* description, Float_t mean, Float_t sigma)
199 {
200   AliTRDCalPad *calPad = new AliTRDCalPad(shortName, description);
201   for (Int_t det=0; det<AliTRDgeometry::kNdet; ++det)
202     {
203       vav[det]=0.0;
204       Int_t nb = 0;
205       AliTRDCalROC *calROC = calPad->GetCalROC(det);
206       for (Int_t channel=0; channel<calROC->GetNchannels(); ++channel){
207         Double_t value = gRandom->Gaus(mean,sigma);
208         vav[det] += value;
209         nb++;
210         //cout << "value: " << value << endl;
211         //if(value < 0) calROC->SetValue(channel, 0);
212         calROC->SetValue(channel, value);
213       }
214       vav[det] /= nb;
215       for (Int_t channel=0; channel<calROC->GetNchannels(); ++channel){
216         Double_t value = calROC->GetValue(channel);
217         //cout << "value: " << value << endl;
218         //if(value < 0) calROC->SetValue(channel, 0);
219         calROC->SetValue(channel, value/vav[det]);
220       }
221     }
222   
223   return calPad;
224 }
225 //____________________________________________________________________________________________________
226 TObject* CreatePadObjectbridge(const char* shortName, const char* description)
227 {
228   AliTRDCalPad *calPad = new AliTRDCalPad(shortName, description);
229   Double_t x[2];
230   for (Int_t det=0; det<AliTRDgeometry::kNdet; ++det)
231     {
232   
233       AliTRDCalROC *calROC = calPad->GetCalROC(det);
234       Int_t rowMax = calROC->GetNrows();
235       Int_t colMax = calROC->GetNcols();
236       gainav[det] = 0.0;   
237
238       // premier calcul
239       for(Int_t row = 0; row < rowMax; ++row){
240         for(Int_t col = 0; col < colMax; ++col){
241           x[0] = row;
242           x[1] = col;
243           if(rowMax == 12) {
244             Double_t value = funcpoly12(&x[0])*funcpoly144(&x[1]);
245             calROC->SetValue(col,row,value);
246             gainav[det] += value;
247           }
248           else {
249             Double_t value = funcpoly16(&x[0])*funcpoly144(&x[1]);
250             calROC->SetValue(col,row,value);
251             gainav[det] += value;
252           }
253         }//col
254       }//row
255       gainav[det] /= (rowMax*colMax);
256
257       //deuxieme calcul
258       for(Int_t row = 0; row < rowMax; ++row){
259         for(Int_t col = 0; col < colMax; ++col){
260           x[0] = row;
261           x[1] = col;
262           //Double_t value = funcpoly12(&x[0])*funcpoly144(&x[1]);
263           Double_t value = calROC->GetValue(col,row);
264           calROC->SetValue(col,row,value/gainav[det]);
265           //gainav[det] += value;
266         }//col
267       }//row
268
269     }//det
270   
271   return calPad;
272 }
273 //__________________________________________________________________________________________
274 TObject* CreatePadObjectbridgev(const char* shortName, const char* description)
275 {
276   AliTRDCalPad *calPad = new AliTRDCalPad(shortName, description);
277   Double_t x[2];
278   for (Int_t det=0; det<AliTRDgeometry::kNdet; ++det)
279     {
280   
281       AliTRDCalROC *calROC = calPad->GetCalROC(det);
282       Int_t rowMax = calROC->GetNrows();
283       Int_t colMax = calROC->GetNcols();
284       vav[det] = 0.0;
285
286       // premiere calcul
287       for(Int_t row = 0; row < rowMax; ++row){
288         for(Int_t col = 0; col < colMax; ++col){
289           x[0] = row;
290           x[1] = col;
291           if(rowMax == 12) {
292             Double_t value = funcpoly12v(&x[0])*funcpoly144v(&x[1]);
293             calROC->SetValue(col,row,value);
294             vav[det] += value;
295           }
296           else {
297             Double_t value = funcpoly16v(&x[0])*funcpoly144v(&x[1]);
298             calROC->SetValue(col,row,value);
299             vav[det] += value;
300           }
301         }//col
302       }//row
303       vav[det] /= (rowMax*colMax);
304
305       //deuxieme calcul
306       for(Int_t row = 0; row < rowMax; ++row){
307         for(Int_t col = 0; col < colMax; ++col){
308           x[0] = row;
309           x[1] = col;
310         
311           //Double_t value = funcpoly12v(&x[0])*funcpoly144v(&x[1]);
312           Double_t value = calROC->GetValue(col,row);
313           calROC->SetValue(col,row,value/vav[det]);
314         
315         }//col
316       }//row
317
318     }//det
319   
320   return calPad;
321 }
322 //___________________________________________________________________________________________________
323 TObject* CreateDetObject(const char* shortName, const char* description, Float_t value)
324 {
325   AliTRDCalDet *object = new AliTRDCalDet(shortName, description);
326   for (Int_t det=0; det<AliTRDgeometry::kNdet; ++det)
327     object->SetValue(det, value);
328   return object;
329 }
330 //___________________________________________________________________________________________________
331 TObject* CreateDetObjectRandom(const char* shortName, const char* description, Float_t mean, Float_t sigma)
332 {
333   AliTRDCalDet *object = new AliTRDCalDet(shortName, description);
334   for (Int_t det=0; det<AliTRDgeometry::kNdet; ++det){
335     Double_t value = gRandom->Gaus(mean,sigma);
336     //cout << "value: " << value << endl;
337     object->SetValue(det, value);
338     //else cout << "value negative!" << endl;
339   }
340   return object;
341 }
342 //___________________________________________________________________________________________________
343 TObject* CreateDetObjectRandomg(const char* shortName, const char* description, Float_t mean, Float_t sigma)
344 {
345   AliTRDCalDet *object = new AliTRDCalDet(shortName, description);
346   for (Int_t det=0; det<AliTRDgeometry::kNdet; ++det){
347     Double_t value = gRandom->Gaus(mean,sigma);
348     //cout << "value: " << value << endl;
349     object->SetValue(det, value/gainav[det]);
350     //else cout << "value negative!" << endl;
351   }
352   return object;
353 }
354 //___________________________________________________________________________________________________
355 TObject* CreateDetObjectRandomv(const char* shortName, const char* description, Float_t mean, Float_t sigma)
356 {
357   AliTRDCalDet *object = new AliTRDCalDet(shortName, description);
358   for (Int_t det=0; det<AliTRDgeometry::kNdet; ++det){
359     Double_t value = gRandom->Gaus(mean,sigma);
360     //cout << "value: " << value << endl;
361     object->SetValue(det, value/vav[det]);
362     //else cout << "value negative!" << endl;
363   }
364   return object;
365 }
366 //___________________________________________________________________________________________________
367 AliCDBMetaData* CreateMetaObject(const char* objectClassName)
368 {
369   AliCDBMetaData *md1= new AliCDBMetaData(); 
370   md1->SetObjectClassName(objectClassName);
371   md1->SetResponsible("Raphaelle Bailhache");
372   //md1->SetBeamPeriod(1);
373   //md1->SetAliRootVersion("head"); //root version
374   md1->SetComment("residual database TRD");
375   
376   return md1;
377 }
378 //___________________________________________________________________________________________________
379 void StoreObject(const char* cdbPath, TObject* object, AliCDBMetaData* metaData)
380 {
381   AliCDBId id1(cdbPath, gkDummyRun, 999999999); 
382   id1.SetVersion(1);
383   gStorLoc->Put(object, id1, metaData); 
384 }
385 //___________________________________________________________________________________________________
386 void AliTRDCreate(Bool_t residual)
387 {
388  
389   //************************random generator******************************************
390   TDatime *datime = new TDatime();
391   Int_t time = datime->GetTime();
392   Int_t date = datime->GetDate();
393   Int_t pid  = gSystem->GetPid();
394   delete datime;
395   Int_t seed = TMath::Abs(10000 * pid + time - date);
396   gRandom->SetSeed(seed); 
397
398   //*************************************************************************
399
400   cout << endl << "TRD :: Creating dummy CDB with event number " << gkDummyRun << endl;
401   
402   AliCDBManager *man = AliCDBManager::Instance();
403   gStorLoc = man->GetStorage("local://.");
404   if (!gStorLoc)
405     return;
406
407   TObject* obj = 0;
408   AliCDBMetaData* metaData = 0;
409
410   if(!residual)
411     {
412
413       //Pad////////////////////////////////////////////////////////////////////
414
415       metaData = CreateMetaObject("AliTRDCalPad");
416       
417       obj = CreatePadObjectT0Random("LocalT0","T0 (local variations)", 0, 0.2);
418       StoreObject("TRD/Calib/LocalT0", obj, metaData);
419
420       obj = CreatePadObjectbridge("LocalGainFactor","GainFactor (local variations)");
421       StoreObject("TRD/Calib/LocalGainFactor", obj, metaData);
422
423       obj = CreatePadObjectbridgev("LocalVdrift","TRD drift velocities (local variations)");
424       StoreObject("TRD/Calib/LocalVdrift", obj, metaData);
425       
426       //Det//////////////////////////////////////////////////////////////////
427       
428       metaData = CreateMetaObject("AliTRDCalDet");
429       
430       obj = CreateDetObjectRandom("ChamberVdrift","TRD drift velocities (detector value)", 1.5, 0.08);
431       StoreObject("TRD/Calib/ChamberVdrift", obj, metaData);
432
433       obj = CreateDetT0Object("ChamberT0","T0 (detector value)");
434       StoreObject("TRD/Calib/ChamberT0", obj, metaData);
435       
436       obj = CreateDetObjectRandom("ChamberGainFactor","GainFactor (detector value)", 1.0, 0.18);
437       StoreObject("TRD/Calib/ChamberGainFactor", obj, metaData);
438
439     }
440   else
441     {
442
443       //Pad////////////////////////////////////////////////////////////////////
444   
445       metaData = CreateMetaObject("AliTRDCalPad");
446       
447       obj = CreatePadObjectRandomv("LocalVdrift","TRD drift velocities (local variations)", 1.5, 0.015);
448       StoreObject("TRD/Calib/LocalVdrift", obj, metaData);
449   
450       obj = CreatePadObjectT0Random("LocalT0","T0 (local variations)", 0, 0.02);
451       StoreObject("TRD/Calib/LocalT0", obj, metaData);
452
453       obj = CreatePadObjectRandom("LocalGainFactor","GainFactor (local variations)", 1, 0.01);
454       StoreObject("TRD/Calib/LocalGainFactor", obj, metaData);
455   
456       //Det//////////////////////////////////////////////////////////////////
457       
458       metaData = CreateMetaObject("AliTRDCalDet");
459   
460       obj = CreateDetT0Object("ChamberT0","T0 (detector value)");
461       StoreObject("TRD/Calib/ChamberT0", obj, metaData);
462       
463     }  
464
465 }