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 }