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
  38 #include "AliESDEvent.h"
  39 #include "AliESDtrack.h"
  40 #include "AliESDfriend.h"
  41
  42 #include "AliTPCcalibTime.h"
  43 #include "AliTPCseed.h"
  44 #include "AliTPCcalibCalib.h"
  45
  46 #include "TObjArray.h"
  47 #include "TString.h"
  48
  49 #include "THnSparse.h"
  50 #include "TGraphErrors.h"
  51
  52 #include <cstdlib>
  53 #include <cerrno>
  54
  55 #include "AliHLTReadoutList.h"
  56
  57 ClassImp(AliHLTTPCCalibTimeComponent) // ROOT macro for the implementation of ROOT specific class methods
  58
  59 AliHLTTPCCalibTimeComponent::AliHLTTPCCalibTimeComponent()
  60   :
  61   fCalibTime(NULL),
  62   fESDevent(NULL),
  63   fESDtrack(NULL),
  64   fESDfriendTrack(NULL),
  65   fSeedArray(NULL),
  66   fMinPartition(5),
  67   fMaxPartition(0),
  68   fMinSlice(35),
  69   fMaxSlice(0),
  70   fSpecification(0) ,
  71   fEnableAnalysis(kTRUE)
  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 AliHLTTPCCalibTimeComponent::~AliHLTTPCCalibTimeComponent() {
  81 // see header file for class documentation
  82 }
  83
  84
  85 const char* AliHLTTPCCalibTimeComponent::GetComponentID() {
  86 // see header file for class documentation
  87
  88   return "TPCCalibTime";
  89 }
  90
  91 void AliHLTTPCCalibTimeComponent::GetInputDataTypes( vector<AliHLTComponentDataType>& list) {
  92 // see header file for class documentation
  93
  94   list.clear();
  95   list.push_back( kAliHLTDataTypeTObjArray|kAliHLTDataOriginTPC ); // output of TPCCalibSeedMaker
  96   list.push_back( kAliHLTDataTypeESDObject|kAliHLTDataOriginOut ); // output of global esd converter
  97 }
  98
  99 AliHLTComponentDataType AliHLTTPCCalibTimeComponent::GetOutputDataType() {
 100 // see header file for class documentation
 101
 102   return AliHLTTPCDefinitions::fgkCalibCEDataType|kAliHLTDataOriginOut;
 103 }
 104
 105 void AliHLTTPCCalibTimeComponent::GetOutputDataSize( unsigned long& constBase, double& inputMultiplier ) {
 106 // see header file for class documentation
 107
 108   constBase = 20000;
 109   inputMultiplier = (2.0); // to be estimated
 110 }
 111
 112 AliHLTComponent* AliHLTTPCCalibTimeComponent::Spawn() {
 113 // see header file for class documentation
 114
 115   return new AliHLTTPCCalibTimeComponent();
 116 }
 117
 118
 119 Int_t AliHLTTPCCalibTimeComponent::ScanArgument( Int_t argc, const char** argv ) {
 120 // see header file for class documentation
 121
 122   Int_t iResult = 0;
 123   TString argument = "";
 124   TString parameter = "";
 125
 126   if(!argc) return -EINVAL;
 127
 128   argument = argv[iResult];
 129
 130   if(argument.IsNull()) return -EINVAL;
 131
 132   if( argument.CompareTo("-enable-analysis") == 0 ){
 133     HLTInfo( "Analysis before shipping data to FXS enabled." );
 134     fEnableAnalysis = kTRUE;
 135   }
 136   else {
 137     iResult = -EINVAL;
 138   }
 139   return iResult;
 140 }
 141
 142 Int_t AliHLTTPCCalibTimeComponent::InitCalibration() {
 143 // see header file for class documentation
 144
 145   if(fCalibTime) return EINPROGRESS;
 146   //fCalibTime = new AliTPCcalibTime();
 147   //fCalibTime = new AliTPCcalibTime("calibTime","time dependent Vdrift calibration",-2, 2, 1);
 148
 149   return 0;
 150 }
 151
 152 Int_t AliHLTTPCCalibTimeComponent::DeinitCalibration() {
 153 // see header file for class documentation
 154
 155   if(fCalibTime) delete fCalibTime; fCalibTime = NULL;
 156
 157   return 0;
 158 }
 159
 160 Int_t AliHLTTPCCalibTimeComponent::ProcessCalibration( const AliHLTComponentEventData& /*evtData*/,  AliHLTComponentTriggerData& /*trigData*/ ){
 161 // see header file for class documentation
 162
 163   if(GetFirstInputBlock( kAliHLTDataTypeSOR ) || GetFirstInputBlock( kAliHLTDataTypeEOR )) return 0;
 164
 165   TObject *iter = NULL;
 166
 167   //--------------- output over TObjArray of AliTPCseed objects (output of TPCSeedMaker) -------------------//
 168
 169   for(iter = (TObject*)GetFirstInputObject(kAliHLTDataTypeTObjArray|kAliHLTDataOriginTPC); iter != NULL; iter = (TObject*)GetNextInputObject()){
 170
 171       if(GetDataType(iter) != (kAliHLTDataTypeTObjArray | kAliHLTDataOriginTPC)) continue;
 172       fSeedArray = dynamic_cast<TObjArray*>(iter);
 173    }
 174
 175
 176   //----------- loop over output of global esd converter ----------------//
 177
 178   for(iter = (TObject*)GetFirstInputObject(kAliHLTDataTypeESDObject | kAliHLTDataOriginOut); iter != NULL; iter = (TObject*)GetNextInputObject()){
 179
 180     if(GetDataType(iter) != (kAliHLTDataTypeESDObject | kAliHLTDataOriginOut)) continue;
 181
 182     AliHLTUInt8_t slice     = AliHLTTPCDefinitions::GetMinSliceNr(GetSpecification(iter));
 183     AliHLTUInt8_t partition = AliHLTTPCDefinitions::GetMinPatchNr(GetSpecification(iter));
 184
 185     if( partition < fMinPartition ) fMinPartition = partition;
 186     if( partition > fMaxPartition ) fMaxPartition = partition;
 187     if( slice < fMinSlice ) fMinSlice = slice;
 188     if( slice > fMaxSlice ) fMaxSlice = slice;
 189
 190     fESDevent = dynamic_cast<AliESDEvent*>(iter);
 191     fESDevent->CreateStdContent();
 192
 193     AliESDfriend *f = new AliESDfriend();
 194     fESDevent->AddObject(f);
 195
 196     //fESDevent->SetTimeStamp(1256910155);
 197     //fESDevent->SetRunNumber(84714);
 198     //HLTInfo("time stamp and event number -----: %d, %d\n", fESDevent->GetTimeStamp(), fESDevent->GetRunNumber());
 199     //fCalibTime->UpdateEventInfo(fESDevent);
 200
 201     HLTDebug("# Seeds: %i\n", fSeedArray->GetEntriesFast());
 202
 203     AliTPCcalibCalib *cal = new AliTPCcalibCalib();
 204
 205     for(Int_t i=0; i<fSeedArray->GetEntriesFast(); i++){
 206
 207         AliTPCseed *seed = (AliTPCseed*)fSeedArray->UncheckedAt(i);
 208         if(!seed) continue;
 209
 210         fESDtrack = fESDevent->GetTrack(i);
 211         if(!fESDtrack) continue;
 212         fESDfriendTrack = const_cast<AliESDfriendTrack*>(fESDtrack->GetFriendTrack());
 213         if(!fESDfriendTrack) continue;
 214
 215         cal->RefitTrack(fESDtrack, seed, GetBz());
 216
 217         AliTPCseed *seedCopy = new AliTPCseed(*seed, kTRUE);
 218         //fESDtrack->AddCalibObject(seedCopy);
 219         fESDfriendTrack->AddCalibObject(seedCopy);
 220
 221     }
 222   }
 223
 224   if(!fCalibTime) {
 225     Int_t startTime=fESDevent->GetTimeStamp()-60*60*1;  //Start time one hour before first event, will make precise cuts later.
 226     Int_t   endTime=fESDevent->GetTimeStamp()+60*60*23; //End time 23 hours after first event.
 227     fCalibTime = new AliTPCcalibTime("calibTime","time dependent Vdrift calibration", startTime, endTime, 20*60);
 228     printf("fCalibTime=%i, startTime=%i, endTime=%i\n", fCalibTime!=0, startTime, endTime);
 229   }
 230
 231   fCalibTime->UpdateEventInfo(fESDevent);
 232   fCalibTime->Process(fESDevent);
 233
 234
 235   fSpecification = AliHLTTPCDefinitions::EncodeDataSpecification( fMinSlice, fMaxSlice, fMinPartition, fMaxPartition );
 236   PushBack( (TObject*)fCalibTime, AliHLTTPCDefinitions::fgkCalibCEDataType | kAliHLTDataOriginOut, fSpecification);
 237
 238   return 0;
 239 }
 240
 241 Int_t AliHLTTPCCalibTimeComponent::ShipDataToFXS( const AliHLTComponentEventData& /*evtData*/,  AliHLTComponentTriggerData& /*trigData*/ ){
 242 // see header file for class documentation
 243
 244   HLTInfo("Shipping data to FXS...\n");
 245
 246   if(fEnableAnalysis) fCalibTime->Analyze();
 247
 248   THnSparse* addHist=fCalibTime->GetHistoDrift("all");
 249   if(!addHist) return -1;
 250
 251   //Identifying used range of histogram
 252   Int_t startTimeBin=0;
 253   Int_t endTimeBin=0;
 254   TH1D* histoTime=addHist->Projection(0);
 255   if(histoTime) {
 256     startTimeBin=histoTime->FindFirstBinAbove(0);
 257     endTimeBin  =histoTime->FindLastBinAbove(0);
 258     printf("startTimeBin=%i endTimeBin=%i\n", startTimeBin, endTimeBin);
 259     printf("startTimeBinCentre=%f endTimeBinCentre=%f\n", histoTime->GetBinCenter(startTimeBin), histoTime->GetBinCenter(endTimeBin));
 260     printf("startTimeBinWidth=%f endTimeBinWidth=%f\n", histoTime->GetBinWidth(startTimeBin), histoTime->GetBinWidth(endTimeBin));
 261     delete histoTime;
 262     histoTime=0;
 263   }
 264
 265   Int_t startPtBin=0;
 266   Int_t endPtBin=0;
 267   TH1D* histoPt=addHist->Projection(1);
 268   if(histoPt) {
 269     startPtBin=histoPt->FindFirstBinAbove(0);
 270     endPtBin  =histoPt->FindLastBinAbove(0);
 271     printf("startPtBin=%i endPtBin=%i\n", startPtBin, endPtBin);
 272     printf("startPtBinCentre=%f endPtBinCentre=%f\n", histoPt->GetBinCenter(startPtBin), histoPt->GetBinCenter(endPtBin));
 273     printf("startPtinWidth=%f endPtBinWidth=%f\n", histoPt->GetBinWidth(startPtBin), histoPt->GetBinWidth(endPtBin));
 274     delete histoPt;
 275     histoPt=0;
 276   }
 277
 278   Int_t startVdBin=0;
 279   Int_t endVdBin=0;
 280   TH1D* histoVd=addHist->Projection(2);
 281   if(histoVd) {
 282     startVdBin=histoVd->FindFirstBinAbove(0);
 283     endVdBin  =histoVd->FindLastBinAbove(0);
 284     printf("startVdBin=%i endVdBin=%i\n", startVdBin, endVdBin);
 285     printf("startVdBinCentre=%f endVdBinCentre=%f\n", histoVd->GetBinCenter(startVdBin), histoVd->GetBinCenter(endVdBin));
 286     printf("startVdBinWidth=%f endVdBinWidth=%f\n", histoVd->GetBinWidth(startVdBin), histoVd->GetBinWidth(endVdBin));
 287     delete histoVd;
 288     histoVd=0;
 289   }
 290
 291   TH1D* histoRun=addHist->Projection(3);
 292   Int_t startRunBin=0;
 293   Int_t endRunBin=0;
 294   if(histoRun) {
 295     startRunBin=histoRun->FindFirstBinAbove(0);
 296     endRunBin  =histoRun->FindLastBinAbove(0);
 297     printf("startRunBin=%i endRunBin=%i\n", startRunBin, endRunBin);
 298     printf("startRunBinCentre=%f endRunBinCentre=%f\n", histoRun->GetBinCenter(startRunBin), histoRun->GetBinCenter(endRunBin));
 299     printf("startRunBinWidth=%f endRunBinWidth=%f\n", histoRun->GetBinWidth(startRunBin), histoRun->GetBinWidth(endRunBin));
 300     delete histoRun;
 301     histoRun=0;
 302   }
 303
 304   TObjArray* vdriftArray = new TObjArray();
 305   if(!vdriftArray) return -2;
 306
 307   TObjArray* array=fCalibTime->GetHistoDrift();
 308   if(!array) return -3;
 309
 310   TIterator* iterator = array->MakeIterator();
 311   if(!iterator) return -4;
 312
 313   iterator->Reset();
 314   THnSparse* hist=NULL;
 315   while((hist=(THnSparseF*)iterator->Next())){
 316     if(!hist) continue;
 317     hist->Print();
 318     hist->GetAxis(0)->SetRange(startTimeBin, endTimeBin);
 319     hist->GetAxis(1)->SetRange(startPtBin,   endPtBin);
 320     hist->GetAxis(0)->SetRange(startVdBin,   endVdBin);
 321     hist->GetAxis(3)->SetRange(startRunBin,  endRunBin);
 322     TString name=hist->GetName();
 323     Int_t dim[4]={0,1,2,3};
 324     THnSparse* newHist=hist->Projection(4,dim);
 325     newHist->SetName(name);
 326     vdriftArray->Add(newHist);
 327     TGraphErrors* graph=AliTPCcalibBase::FitSlices(newHist,2,0,400,100,0.05,0.95, kTRUE);
 328     printf("name=%s graph=%i\n", name.Data(), graph==0);
 329     if(!graph || !graph->GetN()) continue;
 330     printf("name=%s graph=%i, N=%i\n", name.Data(), graph==0, graph->GetN());
 331     Int_t pos=name.Index("_");
 332     name=name(pos,name.Capacity()-pos);
 333     TString graphName=graph->ClassName();
 334     graphName+=name;
 335     graphName.ToUpper();
 336     graph->SetName(graphName);
 337     printf("name=%s\n", graphName.Data());
 338     vdriftArray->Add(graph);
 339
 340     //Currently, AliSplineFits can not be given names...
 341     //AliSplineFit* fit=new AliSplineFit();
 342     //fit->SetGraph(graph);
 343     //fit->SetMinPoints(graph->GetN()+1);
 344     //fit->InitKnots(graph,2,0,0.001);
 345     //fit->SplineFit(0);
 346     //TString fiName=fit->ClassName();
 347     //fiName+=type;
 348     //fiName+=trigger;
 349     //fiName.ToUpper();
 350     //fit->SetName(fiName.Data());
 351     //printf("name=%s\n", fiName.Data());
 352     //vdriftArray->Add(fit);
 353   }
 354   THnSparse* laserHist=NULL;
 355   TGraphErrors* laserGraph=NULL;
 356   TString laserName="";
 357
 358   //Histograms and graphs for A side lasers
 359   laserHist=fCalibTime->GetHistVdriftLaserA(1);
 360   if(laserHist){
 361     laserName=laserHist->ClassName();
 362     laserName+="_MEAN_DRIFT_LASER_ALL_A";
 363     laserName.ToUpper();
 364     laserHist->SetName(laserName);
 365     vdriftArray->Add(laserHist);
 366     laserGraph=AliTPCcalibBase::FitSlices(laserHist,2,0,400,100,0.05,0.95, kTRUE);
 367     if(laserGraph && laserGraph->GetN()){
 368       laserName=laserGraph->GetName();
 369       laserName+="_MEAN_DRIFT_LASER_ALL_A";
 370       laserName.ToUpper();
 371       laserGraph->SetName(laserName);
 372       vdriftArray->Add(laserGraph);
 373     }
 374   }
 375
 376   //Histograms and graphs for C side lasers
 377   laserHist=fCalibTime->GetHistVdriftLaserC(1);
 378   if(laserHist){
 379     laserName=laserHist->ClassName();
 380     laserName+="_MEAN_DRIFT_LASER_ALL_C";
 381     laserName.ToUpper();
 382     laserHist->SetName(laserName);
 383     vdriftArray->Add(laserHist);
 384     laserGraph=AliTPCcalibBase::FitSlices(laserHist,2,0,400,100,0.05,0.95, kTRUE);
 385     if(laserGraph && laserGraph->GetN()){
 386       laserName=laserGraph->GetName();
 387       laserName+="_MEAN_DRIFT_LASER_ALL_C";
 388       laserName.ToUpper();
 389       laserGraph->SetName(laserName);
 390       vdriftArray->Add(laserGraph);
 391     }
 392   }
 393
 394   //Meatdata set in off-line...
 395   //AliCDBMetaData *metaData= new AliCDBMetaData();
 396   //metaData->SetObjectClassName("TObjArray");
 397   //metaData->SetResponsible("Dag Toppe Larsen");
 398   //metaData->SetBeamPeriod(1);
 399   //metaData->SetAliRootVersion("05-25-01"); //root version
 400   //metaData->SetComment("Calibration of the time dependence of the drift velocity due to pressure and temperature changes");
 401   //AliCDBId* id1=NULL;
 402   //if(end) id1=new AliCDBId("TPC/Calib/TimeDrift", runNumber, end);
 403   //else    id1=new AliCDBId("TPC/Calib/TimeDrift", runNumber, runNumber);
 404   //AliCDBStorage* gStorage = AliCDBManager::Instance()->GetStorage("local://$ALICE_ROOT/OCDB");
 405   //gStorage->Put(vdriftArray, (*id1), metaData);
 406   //printf("done runNumber=%i, end=%i\n", runNumber, end);
 407
 408   static AliHLTReadoutList rdList(AliHLTReadoutList::kTPC);
 409   PushToFXS( (TObject*)vdriftArray, "TPC", "Time", rdList.Buffer() );
 410   //PushToFXS( (TObject*)vdriftArray, "TPC", "Time");
 411
 412   //Should array be deleted now?
 413   //  if(vdriftArray){
 414   //      vdriftArray.Clear();
 415   //    delete vdriftArray;
 416   //    vdriftArray=0;
 417   //  }
 418
 419   return 0;
 420 }
 421