PWGPP/TRD/AliTRDrecoTask.cxx

   1 ///////////////////////////////////////////////////////\r
   2 //\r
   3 // Basic class for Performance/Calibration TRD tasks\r
   4 // \r
   5 // It performs generic tasks like :\r
   6 //   - data file manegment\r
   7 //   - reference container management\r
   8 //   - debug container management\r
   9 //   - interaction with AliAnalysisManager\r
  10 //   - Plot functor loop\r
  11 //\r
  12 // Author: Alexandru Bercuci <A.Bercuci@gsi.de>, 10/09/2008\r
  13 //\r
  14 /////////////////////////////////////////////////////////\r
  15 \r
  16 #include "TClass.h"\r
  17 #include "TMethod.h"\r
  18 #include "TMethodCall.h"\r
  19 #include "TMethodArg.h"\r
  20 #include "TFile.h"\r
  21 #include "TChain.h"\r
  22 #include "TList.h"\r
  23 #include "TMap.h"\r
  24 #include "TH1.h"\r
  25 #include "TF1.h"\r
  26 #include "TObjArray.h"\r
  27 #include "TDirectory.h"\r
  28 #include "TTreeStream.h"\r
  29 #include "TBox.h"\r
  30 #include "TLatex.h"\r
  31 #include "TVectorT.h"\r
  32 \r
  33 #include "AliLog.h"\r
  34 #include "AliAnalysisTask.h"\r
  35 #include "AliAnalysisManager.h"\r
  36 #include "AliExternalTrackParam.h"\r
  37 \r
  38 #include "info/AliTRDeventInfo.h"\r
  39 #include "AliTRDrecoTask.h"\r
  40 #include "AliTRDtrackV1.h"\r
  41 #include "AliTRDpidUtil.h"\r
  42 \r
  43 ClassImp(AliTRDrecoTask)\r
  44 \r
  45 TList* AliTRDrecoTask::fgTrendPoint(NULL);\r
  46 TTreeSRedirector* AliTRDrecoTask::fgDebugStream(NULL);\r
  47 //_______________________________________________________\r
  48 AliTRDrecoTask::AliTRDrecoTask()\r
  49   : AliAnalysisTaskSE()\r
  50   ,fNRefFigures(0)\r
  51   ,fDets(NULL)\r
  52   ,fContainer(NULL)\r
  53   ,fEvent(NULL)\r
  54   ,fTracks(NULL)\r
  55   ,fClusters(NULL)\r
  56   ,fkClusters(NULL)\r
  57   ,fkTrack(NULL)\r
  58   ,fkMC(NULL)\r
  59   ,fkESD(NULL)\r
  60   ,fSpecies(-6)\r
  61   ,fPt(-1.)\r
  62   ,fPhi(0.)\r
  63   ,fEta(0.)\r
  64   ,fPlotFuncList(NULL)\r
  65   ,fDetFuncList(NULL)\r
  66   ,fRunTerminate(kFALSE)\r
  67 {\r
  68 // Default constructor\r
  69   snprintf(fNameId, 10, "no name");\r
  70 }\r
  71 \r
  72 //_______________________________________________________\r
  73 AliTRDrecoTask::AliTRDrecoTask(const char *name, const char *title)\r
  74   : AliAnalysisTaskSE(name)\r
  75   ,fNRefFigures(0)\r
  76   ,fDets(NULL)\r
  77   ,fContainer(NULL)\r
  78   ,fEvent(NULL)\r
  79   ,fTracks(NULL)\r
  80   ,fClusters(NULL)\r
  81   ,fkClusters(NULL)\r
  82   ,fkTrack(NULL)\r
  83   ,fkMC(NULL)\r
  84   ,fkESD(NULL)\r
  85   ,fSpecies(-6)\r
  86   ,fPt(-1.)\r
  87   ,fPhi(0.)\r
  88   ,fEta(0.)\r
  89   ,fPlotFuncList(NULL)\r
  90   ,fDetFuncList(NULL)\r
  91   ,fRunTerminate(kFALSE)\r
  92 {\r
  93 // Constructor for all derived performance tasks\r
  94 \r
  95   SetTitle(title);\r
  96   snprintf(fNameId, 10, "no name");\r
  97   DefineInput (1, TObjArray::Class()); // track list\r
  98   DefineInput (2, AliTRDeventInfo::Class()); // event info object\r
  99   DefineInput (3, TObjArray::Class()); // cluster list object\r
 100   DefineOutput(1, TObjArray::Class()); // histogram list\r
 101 }\r
 102 \r
 103 //_______________________________________________________\r
 104 AliTRDrecoTask::~AliTRDrecoTask() \r
 105 {\r
 106 \r
 107   // Generic task destructor\r
 108 \r
 109   AliDebug(2, Form(" Ending task %s[%s]", GetName(), GetTitle()));\r
 110   if(fgDebugStream){ \r
 111     delete fgDebugStream;\r
 112     fgDebugStream = NULL;\r
 113   }\r
 114 \r
 115   if(fPlotFuncList){\r
 116     fPlotFuncList->Delete();\r
 117     delete fPlotFuncList;\r
 118     fPlotFuncList = NULL;\r
 119   }\r
 120   if(fDetFuncList){\r
 121     fDetFuncList->Delete();\r
 122     delete fDetFuncList;\r
 123     fDetFuncList = NULL;\r
 124   }\r
 125   \r
 126   if(fDets){\r
 127     if(fDets->IsOwner()) fDets->Delete();\r
 128     delete fDets;\r
 129     fDets = NULL;\r
 130   }\r
 131 \r
 132   if(fContainer && !(AliAnalysisManager::GetAnalysisManager() && AliAnalysisManager::GetAnalysisManager()->IsProofMode())){\r
 133     if(fContainer->IsOwner()) fContainer->Delete();\r
 134     delete fContainer;\r
 135     fContainer = NULL;\r
 136   }\r
 137 \r
 138   if(fgTrendPoint){\r
 139     TFile::Open("TRD.PerformanceTrend.root", "UPDATE");\r
 140     fgTrendPoint->Write();\r
 141     delete fgTrendPoint;\r
 142     fgTrendPoint=NULL;\r
 143     gFile->Close();\r
 144   }\r
 145 }\r
 146 \r
 147 //_______________________________________________________\r
 148 Int_t AliTRDrecoTask::GetNRefFigures() const  \r
 149 { \r
 150   if(!fNRefFigures) AliWarning("No reference plots available.");\r
 151   return fNRefFigures; \r
 152 } \r
 153 \r
 154 //_______________________________________________________\r
 155 void AliTRDrecoTask::UserCreateOutputObjects()\r
 156 {\r
 157   if(!HasFunctorList()) InitFunctorList();\r
 158   fContainer = Histos();\r
 159   PostData(1, fContainer);\r
 160 }\r
 161 \r
 162 //_______________________________________________________\r
 163 void AliTRDrecoTask::UserExec(Option_t *)\r
 164 {\r
 165 // Loop over Plot functors published by particular tasks\r
 166 \r
 167   fTracks   = dynamic_cast<TObjArray *>(GetInputData(1));\r
 168   fEvent    = dynamic_cast<AliTRDeventInfo *>(GetInputData(2));\r
 169   fClusters = dynamic_cast<TObjArray*>(GetInputData(3));\r
 170 \r
 171   if(!fPlotFuncList){\r
 172     AliWarning("No track functor list defined for the task");\r
 173     return;\r
 174   }\r
 175   if(!fTracks) return;\r
 176   if(!fTracks->GetEntriesFast()) return;\r
 177   else AliDebug(2, Form("Tracks[%d] for %s", fTracks->GetEntriesFast(), GetName()));\r
 178 \r
 179   Int_t itrk(-1);\r
 180   AliTRDtrackInfo *trackInfo(NULL);\r
 181   TIter plotIter(fPlotFuncList);\r
 182   TObjArrayIter trackIter(fTracks);\r
 183   while((trackInfo = dynamic_cast<AliTRDtrackInfo*>(trackIter()))){\r
 184     itrk++; fPt=-1; fEta=0.; fPhi=0.; fSpecies=-6;\r
 185     fkMC    = trackInfo->GetMCinfo();\r
 186     fkESD   = trackInfo->GetESDinfo();\r
 187     if((fkTrack = trackInfo->GetTrack())){\r
 188       // cache properties of the track at TRD entrance\r
 189       // check input track status\r
 190       AliExternalTrackParam *tin(NULL);\r
 191       if(!(tin = fkTrack->GetTrackIn())) AliDebug(2, Form("Missing TRD track[%d] :: entry point.", itrk));\r
 192       else {\r
 193         fPt   = tin->Pt();\r
 194         fEta  = tin->Eta();\r
 195         Double_t xyz[3];\r
 196         if(!tin->GetXYZ(xyz)) AliDebug(2, Form("Failed TRD track[%d] :: global track postion", itrk));\r
 197         else fPhi  = TMath::ATan2(xyz[1], xyz[0]);\r
 198         fSpecies= fkTrack->Charge()*(AliTRDpidUtil::Mass2Pid(fkTrack->GetMass())+1);\r
 199       }\r
 200     } else AliDebug(2, Form("Missing TRD track[%d].", itrk));\r
 201 \r
 202     TMethodCall *plot(NULL);\r
 203     plotIter.Reset();\r
 204     while((plot=dynamic_cast<TMethodCall*>(plotIter()))) plot->Execute(this);\r
 205   }\r
 206   if(!fClusters) return;\r
 207   if(!fDetFuncList){\r
 208     AliDebug(1, "No detector functor list defined for task");\r
 209     return;\r
 210   }\r
 211   TIter detIter(fDetFuncList);\r
 212   for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){\r
 213     if(!(fkClusters = (TObjArray*)fClusters->At(idet))) continue;\r
 214     TMethodCall *det(NULL);\r
 215     detIter.Reset();\r
 216     while((det=dynamic_cast<TMethodCall*>(detIter()))) det->Execute(this);\r
 217   }\r
 218 }\r
 219 \r
 220 //_______________________________________________________\r
 221 Bool_t AliTRDrecoTask::GetRefFigure(Int_t /*ifig*/)\r
 222 {\r
 223   AliWarning("Retrieving reference figures not implemented.");\r
 224   return kFALSE;\r
 225 }\r
 226 \r
 227 //_______________________________________________________\r
 228 Bool_t AliTRDrecoTask::PutTrendValue(const Char_t *name, Double_t val)\r
 229 {\r
 230 // Generic publisher for trend values\r
 231 \r
 232   if(!fgTrendPoint){\r
 233     fgTrendPoint = new TList();\r
 234     fgTrendPoint->SetOwner();\r
 235   }\r
 236   fgTrendPoint->AddLast(new TNamed(Form("%s_%s", GetName(), name), Form("%f", val)));\r
 237   return kTRUE;\r
 238 }\r
 239 \r
 240 //_______________________________________________________\r
 241 void AliTRDrecoTask::InitFunctorList()\r
 242 {\r
 243 // Initialize list of functors\r
 244 \r
 245   TClass *c = this->IsA();\r
 246   if(fPlotFuncList) fPlotFuncList->Clear();\r
 247   if(fDetFuncList) fDetFuncList->Clear();\r
 248 \r
 249   TMethod *m(NULL);\r
 250   TIter methIter(c->GetListOfMethods());\r
 251   while((m=dynamic_cast<TMethod*>(methIter()))){\r
 252     TString name(m->GetName());\r
 253     if(name.BeginsWith("Plot")){\r
 254       if(!fPlotFuncList) fPlotFuncList = new TList();\r
 255       fPlotFuncList->AddLast(new TMethodCall(c, (const char*)name, ""));\r
 256     } else if(name.BeginsWith("Det")){\r
 257       if(!fDetFuncList) fDetFuncList = new TList();\r
 258       fDetFuncList->AddLast(new TMethodCall(c, (const char*)name, ""));\r
 259     }\r
 260   }\r
 261 }\r
 262 \r
 263 //_______________________________________________________\r
 264 Bool_t AliTRDrecoTask::Load(const Char_t *file, const Char_t *dir)\r
 265 {\r
 266 // Generic container loader\r
 267 \r
 268   if(!TFile::Open(file)){\r
 269     AliWarning(Form("Couldn't open file %s.", file));\r
 270     return kFALSE;\r
 271   }\r
 272   if(!gFile->cd(dir)){\r
 273     AliWarning(Form("Couldn't cd to %s in %s.", dir, file));\r
 274     return kFALSE;\r
 275   }\r
 276   TObjArray *o = NULL;\r
 277   if(!(o = (TObjArray*)gDirectory->Get(GetName()))){\r
 278     AliWarning("Missing histogram container.");\r
 279     return kFALSE;\r
 280   }\r
 281   fContainer = (TObjArray*)o->Clone(GetName());\r
 282   gFile->Close();\r
 283   return kTRUE;\r
 284 }\r
 285 \r
 286 //________________________________________________________\r
 287 Bool_t AliTRDrecoTask::LoadDetectorMap(const Char_t *file, const Char_t *dir)\r
 288 {\r
 289 // Load detector map.\r
 290 \r
 291   if(!TFile::Open(file)){\r
 292     AliWarning(Form("Couldn't open file %s.", file));\r
 293     return kFALSE;\r
 294   }\r
 295   if(!gFile->cd(dir)){\r
 296     AliWarning(Form("Couldn't cd to %s in %s.", dir, file));\r
 297     return kFALSE;\r
 298   }\r
 299   TObjArray *info = NULL;\r
 300   if(!(info = (TObjArray*)gDirectory->Get("TRDinfoGen"))){\r
 301     AliWarning("Missing TRDinfoGen container.");\r
 302     return kFALSE;\r
 303   }\r
 304   TObjArray *dets = (TObjArray*)info->FindObject("Chambers");\r
 305   if(!dets){\r
 306     AliWarning("Missing detector map from TRDinfoGen results.");\r
 307     info->ls();\r
 308     return kFALSE;\r
 309   }\r
 310   fDets = (TObjArray*)dets->Clone("Chambers");\r
 311   gFile->Close();\r
 312   return kTRUE;\r
 313 }\r
 314 \r
 315 \r
 316 //________________________________________________________\r
 317 Bool_t AliTRDrecoTask::Save(TObjArray * const results){\r
 318   //\r
 319   // Store the output graphs in a ROOT file\r
 320   // Input TObject array will not be written as Key to the file,\r
 321   // only content itself\r
 322   //\r
 323 \r
 324   TDirectory *cwd = gDirectory;\r
 325   if(!TFile::Open(Form("TRD.Result%s.root", GetName()), "RECREATE")) return kFALSE;\r
 326 \r
 327   TIterator *iter = results->MakeIterator();\r
 328   TObject *inObject = NULL, *outObject = NULL;\r
 329   while((inObject = iter->Next())){\r
 330     outObject = inObject->Clone();\r
 331     outObject->Write(NULL, TObject::kSingleKey);\r
 332   }\r
 333   delete iter;\r
 334   gFile->Close(); delete gFile;\r
 335   cwd->cd(); \r
 336   return kTRUE;\r
 337 }\r
 338 \r
 339 //_______________________________________________________\r
 340 Bool_t AliTRDrecoTask::PostProcess()\r
 341 {\r
 342 // To be implemented by particular tasks\r
 343 \r
 344   AliWarning("Post processing of reference histograms not implemented.");\r
 345   return kTRUE;\r
 346 }\r
 347 \r
 348 //_______________________________________________________\r
 349 void AliTRDrecoTask::MakeDetectorPlot(Int_t ly, const Option_t *opt)\r
 350 {\r
 351 // Draw chamber boundaries in eta/phi plots with misalignments\r
 352 // based on info collected by AliTRDinfoGen\r
 353 \r
 354   if(!fDets){\r
 355     AliWarning("Detector map and status not available.");\r
 356     return;\r
 357   }\r
 358   Float_t xmin(0.), xmax(0.);\r
 359   TBox *gdet = new TBox();\r
 360   gdet->SetLineWidth(kBlack);gdet->SetFillColor(kBlack);\r
 361   Int_t style[] = {0, 3003};\r
 362   for(Int_t idet(0); idet<540; idet++){\r
 363     if(idet%6 != ly) continue;\r
 364     TVectorF *det((TVectorF*)fDets->At(idet));\r
 365     if(!det) continue;\r
 366     Int_t iopt = Int_t((*det)[4]);\r
 367     AliDebug(2, Form("det[%03d] 0[%+4.1f %+4.1f] 1[%+4.1f %+4.1f] opt[%d]", idet, (*det)[0], (*det)[1], (*det)[2], (*det)[3], iopt));\r
 368     if(strcmp(opt, "eta")==0){xmin=(*det)[0]; xmax=(*det)[2];}\r
 369     else if(strcmp(opt, "pad")==0){\r
 370       Int_t stk(AliTRDgeometry::GetStack(idet));\r
 371       xmin=-0.6+16*(4-stk)-(stk<2?4:0); xmax=xmin+(stk==2?12:16)-0.2;\r
 372     } else continue;\r
 373     if(iopt==1){\r
 374       gdet->SetFillStyle(style[1]);gdet->SetFillColor(kBlack);\r
 375       gdet->DrawBox(xmin, (*det)[1], xmax, (*det)[3]);\r
 376     } else {\r
 377       gdet->SetFillStyle(style[0]);\r
 378       gdet->DrawBox(xmin, (*det)[1], xmax, (*det)[3]);\r
 379       if(iopt==2){\r
 380         gdet->SetFillStyle(style[1]);gdet->SetFillColor(kBlue);\r
 381         gdet->DrawBox(xmin, (*det)[1], xmax, 0.5*((*det)[3]+(*det)[1]));\r
 382       } else if(iopt==3){\r
 383         gdet->SetFillStyle(style[1]);gdet->SetFillColor(kRed);\r
 384         gdet->DrawBox(xmin, 0.5*((*det)[3]+(*det)[1]), xmax, (*det)[3]);\r
 385       } else if(iopt!=0) AliError(Form("Wrong chmb. status[%d] for det[%03d]", iopt, idet));\r
 386     }\r
 387   }\r
 388   Float_t dsm = TMath::TwoPi()/AliTRDgeometry::kNsector;\r
 389   xmin=0.;\r
 390   if(strcmp(opt, "pad")==0) xmin=38.;\r
 391   TLatex *sm = new TLatex(); sm->SetTextAlign(22);sm->SetTextColor(0); sm->SetTextFont(32);sm->SetTextSize(0.03);\r
 392   for(Int_t is(0); is<AliTRDgeometry::kNsector; is++) sm->DrawLatex(xmin, -TMath::Pi()+(is+0.5)*dsm, Form("%02d", is>=9?(is-9):(is+9)));\r
 393 }\r
 394 \r
 395 \r
 396 //_______________________________________________________\r
 397 void AliTRDrecoTask::MakeSummary()\r
 398 {\r
 399 // To be implemented by particular tasks\r
 400   AliWarning("Summary not available");\r
 401 }\r
 402 \r
 403 //_______________________________________________________\r
 404 void AliTRDrecoTask::SetDebugLevel(Int_t level)\r
 405 {\r
 406 // Generic debug handler\r
 407 \r
 408   AliAnalysisTaskSE::SetDebugLevel(level);\r
 409   if(DebugLevel()>=1){\r
 410     AliInfo(Form("Debug Level for Task %s set to %d", GetName(), level));\r
 411     TDirectory *savedir = gDirectory;\r
 412     fgDebugStream = new TTreeSRedirector("TRD.DebugPerformance.root");\r
 413     savedir->cd();\r
 414   }\r
 415 }\r
 416 \r
 417 //____________________________________________________________________\r
 418 void AliTRDrecoTask::Terminate(Option_t *)\r
 419 {\r
 420   //\r
 421   // Terminate\r
 422   //\r
 423 \r
 424   if(fgDebugStream){ \r
 425     delete fgDebugStream;\r
 426     fgDebugStream = NULL;\r
 427   }\r
 428   fContainer = dynamic_cast<TObjArray *>(GetOutputData(1));\r
 429   if(fContainer && fRunTerminate){\r
 430     PostProcess();\r
 431     MakeSummary();\r
 432   }\r
 433 }\r
 434 \r
 435 //________________________________________________________\r
 436 void AliTRDrecoTask::Adjust(TF1 *f, TH1 * const h)\r
 437 {\r
 438 // Helper function to avoid duplication of code\r
 439 // Make first guesses on the fit parameters\r
 440 \r
 441   // find the intial parameters of the fit !! (thanks George)\r
 442   Int_t nbinsy = Int_t(.5*h->GetNbinsX());\r
 443   Double_t sum = 0.;\r
 444   for(Int_t jbin=nbinsy-4; jbin<=nbinsy+4; jbin++) sum+=h->GetBinContent(jbin); sum/=9.;\r
 445   f->SetParLimits(0, 0., 3.*sum);\r
 446   f->SetParameter(0, .9*sum);\r
 447 \r
 448   f->SetParLimits(1, -.2, .2);\r
 449   f->SetParameter(1, -0.1);\r
 450 \r
 451   f->SetParLimits(2, 0., 4.e-1);\r
 452   f->SetParameter(2, 2.e-2);\r
 453   if(f->GetNpar() <= 4) return;\r
 454 \r
 455   f->SetParLimits(3, 0., sum);\r
 456   f->SetParameter(3, .1*sum);\r
 457 \r
 458   f->SetParLimits(4, -.3, .3);\r
 459   f->SetParameter(4, 0.);\r
 460 \r
 461   f->SetParLimits(5, 0., 1.e2);\r
 462   f->SetParameter(5, 2.e-1);\r
 463 }\r