HLT/TPCLib/calibration/AliHLTTPCCalibTimeComponent.cxx

   1 // $Id$
   2 /**************************************************************************
   3  * This file is property of and copyright by the ALICE HLT Project        *
   4  * ALICE Experiment at CERN, All rights reserved.                         *
   5  *                                                                        *
   6  * Primary Authors: Jochen Thaeder <thaeder@kip.uni-heidelberg.de>        *
   7  *                  for The ALICE HLT Project.                            *
   8  *                                                                        *
   9  * Permission to use, copy, modify and distribute this software and its   *
  10  * documentation strictly for non-commercial purposes is hereby granted   *
  11  * without fee, provided that the above copyright notice appears in all   *
  12  * copies and that both the copyright notice and this permission notice   *
  13  * appear in the supporting documentation. The authors make no claims     *
  14  * about the suitability of this software for any purpose. It is          *
  15  * provided "as is" without express or implied warranty.                  *
  16  **************************************************************************/
  17
  18 /** @file   AliHLTTPCCalibTimeComponent.cxx
  19     @author Kalliopi Kanaki
  20     @date   2009-07-08
  21     @brief  A calibration component for interfacing the offline calculation of TPC drift velocity correction
  22 */
  23
  24 // see header file for class documentation
  25 // or
  26 // refer to README to build package
  27 // or
  28 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  29
  30 #if __GNUC__>= 3
  31 using namespace std;
  32 #endif
  33
  34
  35 #include "AliHLTTPCCalibTimeComponent.h"
  36 #include "AliHLTTPCDefinitions.h"
  37 #include "AliHLTMisc.h"
  38
  39 #include "AliESDEvent.h"
  40 #include "AliESDtrack.h"
  41 #include "AliESDfriend.h"
  42
  43 #include "AliTPCcalibTime.h"
  44 #include "AliTPCcalibCalib.h"
  45 #include "AliTPCseed.h"
  46 #include "AliTPCcalibDB.h"
  47 #include "AliTPCClusterParam.h"
  48
  49 #include "TObjArray.h"
  50 #include "TString.h"
  51 #include "TFile.h"
  52
  53 #include "THnSparse.h"
  54 #include "TGraphErrors.h"
  55
  56 #include <cstdlib>
  57 #include <cerrno>
  58
  59 #include "AliHLTReadoutList.h"
  60
  61 ClassImp(AliHLTTPCCalibTimeComponent) // ROOT macro for the implementation of ROOT specific class methods
  62
  63 AliHLTTPCCalibTimeComponent::AliHLTTPCCalibTimeComponent()
  64   :
  65    fCalibTime(NULL)
  66   ,fCal(NULL)
  67   ,fESDevent(NULL)
  68   ,fESDtrack(NULL)
  69   ,fESDfriend(NULL)
  70   ,fSeedArray(NULL)
  71   ,fOutputSize(50000)
  72 {
  73   // see header file for class documentation
  74   // or
  75   // refer to README to build package
  76   // or
  77   // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
  78 }
  79
  80 const char* AliHLTTPCCalibTimeComponent::fgkOCDBEntry="HLT/ConfigTPC/TPCCalibTime";
  81
  82 AliHLTTPCCalibTimeComponent::~AliHLTTPCCalibTimeComponent() {
  83 // see header file for class documentation
  84 }
  85
  86
  87 const char* AliHLTTPCCalibTimeComponent::GetComponentID() {
  88 // see header file for class documentation
  89
  90   return "TPCCalibTime";
  91 }
  92
  93 void AliHLTTPCCalibTimeComponent::GetInputDataTypes( vector<AliHLTComponentDataType>& list) {
  94 // see header file for class documentation
  95
  96   list.clear();
  97   list.push_back( kAliHLTDataTypeTObjArray|kAliHLTDataOriginTPC ); // output of TPCCalibSeedMaker
  98   list.push_back( kAliHLTDataTypeESDObject|kAliHLTDataOriginOut ); // output of global esd converter
  99 }
 100
 101 AliHLTComponentDataType AliHLTTPCCalibTimeComponent::GetOutputDataType() {
 102 // see header file for class documentation
 103
 104   return AliHLTTPCDefinitions::fgkCalibCEDataType|kAliHLTDataOriginOut;
 105 }
 106
 107 void AliHLTTPCCalibTimeComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier ) {
 108 // see header file for class documentation
 109
 110   constBase = fOutputSize;
 111   inputMultiplier = 0; // to be estimated
 112 }
 113
 114 AliHLTComponent* AliHLTTPCCalibTimeComponent::Spawn() {
 115 // see header file for class documentation
 116
 117   return new AliHLTTPCCalibTimeComponent();
 118 }
 119
 120
 121 Int_t AliHLTTPCCalibTimeComponent::ScanConfigurationArgument( Int_t argc, const char** argv ) {
 122 // see header file for class documentation
 123
 124   if (argc<=0) return 0;
 125   int i=0;
 126   TString argument=argv[i];
 127
 128   // -output-size
 129   if (argument.CompareTo("-output-size")==0) {
 130     if (++i>=argc) return -EPROTO;
 131     argument=argv[i];
 132     fOutputSize=argument.Atof();
 133     return 2;
 134   }
 135   return -EINVAL;
 136 }
 137
 138 Int_t AliHLTTPCCalibTimeComponent::InitCalibration() {
 139 // see header file for class documentation
 140
 141   //AliTPCcalibDB::Instance()->SetRun(84714);
 142   AliTPCcalibDB::Instance()->SetRun(AliHLTMisc::Instance().GetCDBRunNo());
 143   AliTPCcalibDB::Instance()->GetClusterParam()->SetInstance(AliTPCcalibDB::Instance()->GetClusterParam());
 144
 145
 146 //   AliTPCcalibDB *calib = AliTPCcalibDB::Instance();
 147 //
 148 //   if(!calib){
 149 //     HLTError("AliTPCcalibDB does not exist");
 150 //     return -ENOENT;
 151 //   }
 152 //
 153 //   AliTPCClusterParam *clusPar = calib->GetClusterParam();
 154 //   if(!clusPar){
 155 //     HLTError("OCDB entry TPC/Calib/ClusterParam (AliTPCcalibDB::GetClusterParam()) is not available.");
 156 //     return -ENOENT;
 157 //   }
 158
 159   // first configure the default
 160   int iResult=0;
 161   if (iResult>=0) iResult=ConfigureFromCDBTObjString(fgkOCDBEntry);
 162
 163   // configure from the command line parameters if specified
 164   //if (iResult>=0 && argc>0)  iResult=ConfigureFromArgumentString(argc, argv);
 165
 166   if(fCalibTime) return EINPROGRESS;
 167   fCal = new AliTPCcalibCalib();
 168
 169   return iResult;
 170
 171 }
 172
 173 Int_t AliHLTTPCCalibTimeComponent::DeinitCalibration() {
 174 // see header file for class documentation
 175
 176   if(fCalibTime) delete fCalibTime; fCalibTime = NULL;
 177   if(fCal)       delete fCal;       fCal       = NULL;
 178   //if(fESDfriend) delete fESDfriend; fESDfriend = NULL;
 179
 180   return 0;
 181 }
 182
 183 int AliHLTTPCCalibTimeComponent::Reconfigure(const char* cdbEntry, const char* /*chainId*/){
 184 // see header file for class documentation
 185
 186   // configure from the specified antry or the default one
 187   const char* entry=cdbEntry;
 188   if (!entry || entry[0]==0) {
 189      entry=fgkOCDBEntry;
 190   }
 191
 192   return ConfigureFromCDBTObjString(entry);
 193 }
 194
 195 Int_t AliHLTTPCCalibTimeComponent::ProcessCalibration( const AliHLTComponentEventData& /*evtData*/,  AliHLTComponentTriggerData& /*trigData*/ ){
 196 // see header file for class documentation
 197
 198   if(GetFirstInputBlock( kAliHLTDataTypeSOR ) || GetFirstInputBlock( kAliHLTDataTypeEOR )) return 0;
 199
 200   TObject *iter = NULL;
 201
 202   //--------------- output over TObjArray of AliTPCseed objects (output of TPCSeedMaker) -------------------//
 203
 204   // A previous component in the chain (TPCSeedMaker) has processed the TPC clusters and tracks and created a TObjArray of AliTPCseed objects
 205   // In this loop the iterator accesses this array stored in memory, in order to use it in the next loop over the AliESDevent of the HLT
 206
 207   for(iter = (TObject*)GetFirstInputObject(kAliHLTDataTypeTObjArray|kAliHLTDataOriginTPC); iter != NULL; iter = (TObject*)GetNextInputObject()){
 208
 209       if(GetDataType(iter) != (kAliHLTDataTypeTObjArray | kAliHLTDataOriginTPC)) continue;
 210       fSeedArray = dynamic_cast<TObjArray*>(iter);
 211    }
 212
 213
 214   //----------- loop over output of global esd converter ----------------//
 215
 216   // In this loop we access the AliESDevent that was produced by the HLT and is stored in memory. There should exist 1 object
 217   // of type kAliHLTDataTypeESDObject per event.
 218
 219   for(iter = (TObject*)GetFirstInputObject(kAliHLTDataTypeESDObject | kAliHLTDataOriginOut); iter != NULL; iter = (TObject*)GetNextInputObject()){
 220
 221     if(GetDataType(iter) != (kAliHLTDataTypeESDObject | kAliHLTDataOriginOut)) continue;
 222
 223     fESDevent = dynamic_cast<AliESDEvent*>(iter);
 224     fESDevent->GetStdContent();
 225
 226     //fESDevent->SetTimeStamp(1256910155);
 227     //fESDevent->SetRunNumber(0);
 228     //fESDevent->SetRunNumber(84714);
 229
 230     HLTDebug("# Seeds: %i\n", fSeedArray->GetEntriesFast()); // access of the info from the previous loop over the AliTPCseed array
 231
 232     fCal->UpdateEventInfo(fESDevent);
 233
 234     for(Int_t i=0; i<fSeedArray->GetEntriesFast(); i++){  // loop over TObjArray
 235
 236         AliTPCseed *seed = (AliTPCseed*)fSeedArray->UncheckedAt(i);
 237         if(!seed) continue;
 238
 239         fESDtrack = fESDevent->GetTrack(i);
 240         if(!fESDtrack) continue;
 241
 242         fCal->RefitTrack(fESDtrack, seed, GetBz()); // update AliESDtrack and AliTPCseed info, acccording to Marian's request
 243
 244         AliTPCseed *seedCopy = new AliTPCseed(*seed, kTRUE);
 245         fESDtrack->AddCalibObject(seedCopy);  // add the AliTPCseed as a friend track to the AliESDtrack (to be accessed in TPC/AliTPCcalibTime.cxx)
 246
 247         //fESDfriendTrack = const_cast<AliESDfriendTrack*>(fESDtrack->GetFriendTrack());
 248    }
 249   }
 250
 251   if(!fCalibTime){ // create the calibration object that will call the offline functions
 252
 253      Int_t startTime = fESDevent->GetTimeStamp()-60*60*1;  //Start time one hour before first event, will make precise cuts later.
 254      Int_t   endTime = fESDevent->GetTimeStamp()+60*60*23; //End time 23 hours after first event.
 255      fCalibTime = new AliTPCcalibTime("calibTime","time dependent Vdrift calibration", startTime, endTime, 20*60);
 256      fCalibTime->SetStreamLevel(20);
 257      fCalibTime->SetDebugLevel(20);
 258      printf("fCalibTime = %i, startTime = %i, endTime = %i \n", fCalibTime!=0, startTime, endTime);
 259   }
 260
 261   fESDfriend = new AliESDfriend();
 262   fESDevent->GetESDfriend(fESDfriend);
 263   fESDevent->SetESDfriend(fESDfriend);
 264   fESDevent->AddObject(fESDfriend);
 265   // create the AliESDfriend and add it to the event, now both the friend tracks and the friends are available for the offline functions to be called
 266
 267   fCalibTime->UpdateEventInfo(fESDevent); // needed for getting the run number and time stamp information correct on the offline side
 268   fCalibTime->Process(fESDevent);         // first offline function called
 269
 270   // delete fESDfriend;
 271
 272   //PushBack( (TObject*)fCalibTime, AliHLTTPCDefinitions::fgkCalibCEDataType | kAliHLTDataOriginOut, 0x0);
 273
 274   return 0;
 275 }
 276
 277 Int_t AliHLTTPCCalibTimeComponent::ShipDataToFXS( const AliHLTComponentEventData& /*evtData*/,  AliHLTComponentTriggerData& /*trigData*/ ){
 278 // see header file for class documentation
 279
 280   HLTInfo("Shipping data to FXS...\n");
 281
 282   fCalibTime->Analyze(); // called at the end of the run or event modulo
 283
 284   // the rest of the histogram and graph declarations were copied by Dag as a first attempt to get the start/end time bin "automatically". Perhaps we need some more
 285   // thinking here to avoid copying all these lines that might chain in offline without HLT realizing.
 286
 287   THnSparse* addHist = fCalibTime->GetHistoDrift("all");
 288   if(!addHist) return -1;
 289
 290   //Identifying used range of histogram
 291
 292   Int_t startTimeBin = 0;
 293   Int_t endTimeBin   = 0;
 294
 295   TH1D *histoTime = addHist->Projection(0);
 296   if(histoTime){
 297     startTimeBin = histoTime->FindFirstBinAbove(0);
 298     endTimeBin   = histoTime->FindLastBinAbove(0);
 299     printf("startTimeBin       = %i endTimeBin       = %i\n", startTimeBin, endTimeBin);
 300     printf("startTimeBinCentre = %f endTimeBinCentre = %f\n", histoTime->GetBinCenter(startTimeBin), histoTime->GetBinCenter(endTimeBin));
 301     printf("startTimeBinWidth  = %f endTimeBinWidth  = %f\n", histoTime->GetBinWidth(startTimeBin),  histoTime->GetBinWidth(endTimeBin));
 302     delete histoTime; histoTime = 0;
 303   }
 304
 305   Int_t startPtBin = 0;
 306   Int_t endPtBin   = 0;
 307   TH1D *histoPt = addHist->Projection(1);
 308   if(histoPt){
 309     startPtBin = histoPt->FindFirstBinAbove(0);
 310     endPtBin   = histoPt->FindLastBinAbove(0);
 311     printf("startPtBin       = %i endPtBin       = %i\n", startPtBin, endPtBin);
 312     printf("startPtBinCentre = %f endPtBinCentre = %f\n", histoPt->GetBinCenter(startPtBin), histoPt->GetBinCenter(endPtBin));
 313     printf("startPtinWidth   = %f endPtBinWidth  = %f\n", histoPt->GetBinWidth(startPtBin),  histoPt->GetBinWidth(endPtBin));
 314     delete histoPt; histoPt = 0;
 315   }
 316
 317   Int_t startVdBin = 0;
 318   Int_t endVdBin   = 0;
 319   TH1D *histoVd = addHist->Projection(2);
 320   if(histoVd){
 321     startVdBin = histoVd->FindFirstBinAbove(0);
 322     endVdBin   = histoVd->FindLastBinAbove(0);
 323     printf("startVdBin       = %i endVdBin       = %i\n", startVdBin, endVdBin);
 324     printf("startVdBinCentre = %f endVdBinCentre = %f\n", histoVd->GetBinCenter(startVdBin), histoVd->GetBinCenter(endVdBin));
 325     printf("startVdBinWidth  = %f endVdBinWidth  = %f\n", histoVd->GetBinWidth(startVdBin),  histoVd->GetBinWidth(endVdBin));
 326     delete histoVd; histoVd = 0;
 327   }
 328
 329   Int_t startRunBin = 0;
 330   Int_t endRunBin   = 0;
 331   TH1D *histoRun = addHist->Projection(3);
 332   if(histoRun){
 333     startRunBin = histoRun->FindFirstBinAbove(0);
 334     endRunBin   = histoRun->FindLastBinAbove(0);
 335     printf("startRunBin       = %i endRunBin       = %i\n", startRunBin, endRunBin);
 336     printf("startRunBinCentre = %f endRunBinCentre = %f\n", histoRun->GetBinCenter(startRunBin), histoRun->GetBinCenter(endRunBin));
 337     printf("startRunBinWidth  = %f endRunBinWidth  = %f\n", histoRun->GetBinWidth(startRunBin),  histoRun->GetBinWidth(endRunBin));
 338     delete histoRun; histoRun = 0;
 339   }
 340
 341   TObjArray *vdriftArray = new TObjArray();
 342   if(!vdriftArray) return -2;
 343
 344   TObjArray *array = fCalibTime->GetHistoDrift();
 345   if(!array) return -3;
 346
 347   TIterator *iterator = array->MakeIterator();
 348   if(!iterator) return -4;
 349
 350   iterator->Reset();
 351   THnSparse *hist = NULL;
 352   while((hist = (THnSparseF*)iterator->Next())){
 353
 354          if(!hist) continue;
 355          hist->Print();
 356          hist->GetAxis(0)->SetRange(startTimeBin, endTimeBin);
 357          hist->GetAxis(1)->SetRange(startPtBin,   endPtBin);
 358          hist->GetAxis(0)->SetRange(startVdBin,   endVdBin);
 359          hist->GetAxis(3)->SetRange(startRunBin,  endRunBin);
 360
 361          TString name = hist->GetName();
 362          Int_t dim[4] = {0,1,2,3};
 363          THnSparse *newHist = hist->Projection(4,dim);
 364          newHist->SetName(name);
 365          vdriftArray->Add(newHist);
 366
 367          TGraphErrors *graph = AliTPCcalibBase::FitSlices(newHist,2,0,400,100,0.05,0.95, kTRUE);
 368          printf("name = %s graph = %i\n", name.Data(), graph==0);
 369          if(!graph || !graph->GetN()) continue;
 370          printf("name = %s graph = %i, N = %i\n", name.Data(), graph==0, graph->GetN());
 371          Int_t pos = name.Index("_");
 372          name = name(pos,name.Capacity()-pos);
 373          TString graphName = graph->ClassName();
 374          graphName+=name;
 375          graphName.ToUpper();
 376          graph->SetName(graphName);
 377          printf("name = %s\n", graphName.Data());
 378          vdriftArray->Add(graph);
 379
 380          //Currently, AliSplineFits can not be given names...
 381          //AliSplineFit* fit=new AliSplineFit();
 382          //fit->SetGraph(graph);
 383          //fit->SetMinPoints(graph->GetN()+1);
 384          //fit->InitKnots(graph,2,0,0.001);
 385          //fit->SplineFit(0);
 386          //TString fiName=fit->ClassName();
 387          //fiName+=type;
 388          //fiName+=trigger;
 389          //fiName.ToUpper();
 390          //fit->SetName(fiName.Data());
 391          //printf("name=%s\n", fiName.Data());
 392          //vdriftArray->Add(fit);
 393   }
 394
 395   THnSparse    *laserHist  = NULL;
 396   TGraphErrors *laserGraph = NULL;
 397   TString laserName = "";
 398
 399   //Histograms and graphs for A side lasers
 400   laserHist = fCalibTime->GetHistVdriftLaserA(1);
 401   if(laserHist){
 402
 403      laserName=laserHist->ClassName();
 404      laserName+="_MEAN_DRIFT_LASER_ALL_A";
 405      laserName.ToUpper();
 406      laserHist->SetName(laserName);
 407      vdriftArray->Add(laserHist);
 408      laserGraph=AliTPCcalibBase::FitSlices(laserHist,2,0,400,100,0.05,0.95, kTRUE);
 409      if(laserGraph && laserGraph->GetN()){
 410        laserName=laserGraph->GetName();
 411        laserName+="_MEAN_DRIFT_LASER_ALL_A";
 412        laserName.ToUpper();
 413        laserGraph->SetName(laserName);
 414        vdriftArray->Add(laserGraph);
 415      }
 416   }
 417
 418   //Histograms and graphs for C side lasers
 419   laserHist=fCalibTime->GetHistVdriftLaserC(1);
 420   if(laserHist){
 421      laserName=laserHist->ClassName();
 422      laserName+="_MEAN_DRIFT_LASER_ALL_C";
 423      laserName.ToUpper();
 424      laserHist->SetName(laserName);
 425      vdriftArray->Add(laserHist);
 426      laserGraph=AliTPCcalibBase::FitSlices(laserHist,2,0,400,100,0.05,0.95, kTRUE);
 427      if(laserGraph && laserGraph->GetN()){
 428        laserName=laserGraph->GetName();
 429        laserName+="_MEAN_DRIFT_LASER_ALL_C";
 430        laserName.ToUpper();
 431        laserGraph->SetName(laserName);
 432        vdriftArray->Add(laserGraph);
 433      }
 434   }
 435
 436   //Meatdata set in off-line...
 437   //AliCDBMetaData *metaData= new AliCDBMetaData();
 438   //metaData->SetObjectClassName("TObjArray");
 439   //metaData->SetResponsible("Dag Toppe Larsen");
 440   //metaData->SetBeamPeriod(1);
 441   //metaData->SetAliRootVersion("05-25-01"); //root version
 442   //metaData->SetComment("Calibration of the time dependence of the drift velocity due to pressure and temperature changes");
 443   //AliCDBId* id1=NULL;
 444   //if(end) id1=new AliCDBId("TPC/Calib/TimeDrift", runNumber, end);
 445   //else    id1=new AliCDBId("TPC/Calib/TimeDrift", runNumber, runNumber);
 446   //AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage("local://$ALICE_ROOT/OCDB");
 447   //gStorage->Put(vdriftArray, (*id1), metaData);
 448   //printf("done runNumber=%i, end=%i\n", runNumber, end);
 449
 450   static AliHLTReadoutList rdList(AliHLTReadoutList::kTPC);
 451
 452   // the vdriftArray is pushed to the HLT-FXSsubscriber
 453   PushToFXS( (TObject*)vdriftArray, "TPC", "TIMEDRIFT", rdList.Buffer() );
 454
 455   //PushToFXS( (TObject*)vdriftArray, "TPC", "Time");
 456
 457   TFile *file = TFile::Open("vdrift.root", "RECREATE");
 458   vdriftArray->Write();
 459   file->Close();
 460   delete file;
 461
 462   file = TFile::Open("calibTime.root", "RECREATE");
 463   fCalibTime->Write();
 464   file->Close();
 465   delete file;
 466
 467   //Should array be deleted now?
 468   //  if(vdriftArray){
 469   //      vdriftArray.Clear();
 470   //    delete vdriftArray;
 471   //    vdriftArray=0;
 472   //  }
 473
 474   return 0;
 475 }
 476