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 "TH2.h"\r
  26 #include "TH3.h"\r
  27 #include "TF1.h"\r
  28 #include "TObjArray.h"\r
  29 #include "TDirectory.h"\r
  30 #include "TTreeStream.h"\r
  31 #include "TBox.h"\r
  32 #include "TLatex.h"\r
  33 #include "TVectorT.h"\r
  34 \r
  35 #include "AliLog.h"\r
  36 #include "AliAnalysisTask.h"\r
  37 #include "AliAnalysisManager.h"\r
  38 #include "AliExternalTrackParam.h"\r
  39 \r
  40 #include "info/AliTRDeventInfo.h"\r
  41 #include "info/AliTRDtrendingManager.h"\r
  42 #include "AliTRDrecoTask.h"\r
  43 #include "AliTRDtrackV1.h"\r
  44 #include "AliTRDpidUtil.h"\r
  45 \r
  46 ClassImp(AliTRDrecoTask)\r
  47 \r
  48 Float_t AliTRDrecoTask::fgPt0[AliTRDrecoTask::fgNPt0] = {0.5, 0.8, 1.5, 5};\r
  49 TTreeSRedirector* AliTRDrecoTask::fgDebugStream(NULL);\r
  50 //_______________________________________________________\r
  51 AliTRDrecoTask::AliTRDrecoTask()\r
  52   : AliAnalysisTaskSE()\r
  53   ,fNRefFigures(0)\r
  54   ,fDets(NULL)\r
  55   ,fDetsV(NULL)\r
  56   ,fContainer(NULL)\r
  57   ,fEvent(NULL)\r
  58   ,fTracks(NULL)\r
  59   ,fClusters(NULL)\r
  60   ,fkClusters(NULL)\r
  61   ,fkTrack(NULL)\r
  62   ,fkMC(NULL)\r
  63   ,fkESD(NULL)\r
  64   ,fSpecies(-6)\r
  65   ,fPt(-1.)\r
  66   ,fPhi(0.)\r
  67   ,fEta(0.)\r
  68   ,fPlotFuncList(NULL)\r
  69   ,fDetFuncList(NULL)\r
  70   ,fRunTerminate(kFALSE)\r
  71 {\r
  72 // Default constructor\r
  73   snprintf(fNameId, 10, "no name");\r
  74 }\r
  75 \r
  76 //_______________________________________________________\r
  77 AliTRDrecoTask::AliTRDrecoTask(const char *name, const char *title)\r
  78   : AliAnalysisTaskSE(name)\r
  79   ,fNRefFigures(0)\r
  80   ,fDets(NULL)\r
  81   ,fDetsV(NULL)\r
  82   ,fContainer(NULL)\r
  83   ,fEvent(NULL)\r
  84   ,fTracks(NULL)\r
  85   ,fClusters(NULL)\r
  86   ,fkClusters(NULL)\r
  87   ,fkTrack(NULL)\r
  88   ,fkMC(NULL)\r
  89   ,fkESD(NULL)\r
  90   ,fSpecies(-6)\r
  91   ,fPt(-1.)\r
  92   ,fPhi(0.)\r
  93   ,fEta(0.)\r
  94   ,fPlotFuncList(NULL)\r
  95   ,fDetFuncList(NULL)\r
  96   ,fRunTerminate(kFALSE)\r
  97 {\r
  98 // Constructor for all derived performance tasks\r
  99 \r
 100   SetTitle(title);\r
 101   snprintf(fNameId, 10, "no name");\r
 102   DefineInput (1, TObjArray::Class()); // track list\r
 103   DefineInput (2, AliTRDeventInfo::Class()); // event info object\r
 104   DefineInput (3, TObjArray::Class()); // cluster list object\r
 105   DefineOutput(1, TObjArray::Class()); // histogram list\r
 106 }\r
 107 \r
 108 //_______________________________________________________\r
 109 AliTRDrecoTask::~AliTRDrecoTask() \r
 110 {\r
 111 \r
 112   // Generic task destructor\r
 113 \r
 114   AliDebug(2, Form(" Ending task %s[%s]", GetName(), GetTitle()));\r
 115   if(fgDebugStream){ \r
 116     delete fgDebugStream;\r
 117     fgDebugStream = NULL;\r
 118   }\r
 119 \r
 120   if(fPlotFuncList){\r
 121     fPlotFuncList->Delete();\r
 122     delete fPlotFuncList;\r
 123     fPlotFuncList = NULL;\r
 124   }\r
 125   if(fDetFuncList){\r
 126     fDetFuncList->Delete();\r
 127     delete fDetFuncList;\r
 128     fDetFuncList = NULL;\r
 129   }\r
 130   \r
 131   if(fDets){\r
 132     if(fDets->IsOwner()) fDets->Delete();\r
 133     delete fDets;\r
 134     fDets = NULL;\r
 135   }\r
 136   if(fDetsV) delete fDetsV; fDetsV=NULL;\r
 137 \r
 138   if(fContainer && !(AliAnalysisManager::GetAnalysisManager() && AliAnalysisManager::GetAnalysisManager()->IsProofMode())){\r
 139     if(fContainer->IsOwner()) fContainer->Delete();\r
 140     delete fContainer;\r
 141     fContainer = NULL;\r
 142   }\r
 143 \r
 144 /*  if(fgTrendPoint){\r
 145     TFile::Open("TRD.PerformanceTrend.root", "UPDATE");\r
 146     fgTrendPoint->Write();\r
 147     delete fgTrendPoint;\r
 148     fgTrendPoint=NULL;\r
 149     gFile->Close();\r
 150   }*/\r
 151 }\r
 152 \r
 153 //_______________________________________________________\r
 154 Int_t AliTRDrecoTask::GetNRefFigures() const  \r
 155 { \r
 156   if(!fNRefFigures) AliWarning("No reference plots available.");\r
 157   return fNRefFigures; \r
 158 } \r
 159 \r
 160 //____________________________________________________________________\r
 161 Int_t AliTRDrecoTask::GetPtBinSignificant(Float_t pt)\r
 162 {\r
 163 // Get significant (very low, low, medium, high, very high) pt bin\r
 164 \r
 165   Int_t ipt(0);\r
 166   while(ipt<fgNPt0){\r
 167     if(pt<fgPt0[ipt]) break;\r
 168     ipt++;\r
 169   }\r
 170   return ipt-1;\r
 171 }\r
 172 \r
 173 //_______________________________________________________\r
 174 void AliTRDrecoTask::UserCreateOutputObjects()\r
 175 {\r
 176   if(!HasFunctorList()) InitFunctorList();\r
 177   fContainer = Histos();\r
 178   PostData(1, fContainer);\r
 179 }\r
 180 \r
 181 //_______________________________________________________\r
 182 void AliTRDrecoTask::UserExec(Option_t *)\r
 183 {\r
 184 // Loop over Plot functors published by particular tasks\r
 185 \r
 186   fTracks   = dynamic_cast<TObjArray *>(GetInputData(1));\r
 187   fEvent    = dynamic_cast<AliTRDeventInfo *>(GetInputData(2));\r
 188   fClusters = dynamic_cast<TObjArray*>(GetInputData(3));\r
 189 \r
 190   if(!fPlotFuncList){\r
 191     AliWarning("No track functor list defined for the task");\r
 192     return;\r
 193   }\r
 194   if(!fTracks) return;\r
 195   if(!fTracks->GetEntriesFast()) return;\r
 196   else AliDebug(2, Form("Tracks[%d] for %s", fTracks->GetEntriesFast(), GetName()));\r
 197 \r
 198   Int_t itrk(-1);\r
 199   AliTRDtrackInfo *trackInfo(NULL);\r
 200   TIter plotIter(fPlotFuncList);\r
 201   TObjArrayIter trackIter(fTracks);\r
 202   while((trackInfo = dynamic_cast<AliTRDtrackInfo*>(trackIter()))){\r
 203     itrk++; fPt=-1; fEta=0.; fPhi=0.; fSpecies=-6;\r
 204     fkMC    = trackInfo->GetMCinfo();\r
 205     fkESD   = trackInfo->GetESDinfo();\r
 206     if((fkTrack = trackInfo->GetTrack())){\r
 207       // cache properties of the track at TRD entrance\r
 208       // check input track status\r
 209       AliExternalTrackParam *tin(NULL);\r
 210       if(!(tin = fkTrack->GetTrackIn())) AliDebug(2, Form("Missing TRD track[%d] :: entry point.", itrk));\r
 211       else {\r
 212         fPt   = tin->Pt();\r
 213         fEta  = tin->Eta();\r
 214         Double_t xyz[3];\r
 215         if(!tin->GetXYZ(xyz)) AliDebug(2, Form("Failed TRD track[%d] :: global track postion", itrk));\r
 216         else fPhi  = TMath::ATan2(xyz[1], xyz[0]);\r
 217         fSpecies= fkTrack->Charge()*(AliTRDpidUtil::Mass2Pid(fkTrack->GetMass())+1);\r
 218       }\r
 219     } else AliDebug(2, Form("Missing TRD track[%d].", itrk));\r
 220 \r
 221     TMethodCall *plot(NULL);\r
 222     plotIter.Reset();\r
 223     while((plot=dynamic_cast<TMethodCall*>(plotIter()))) plot->Execute(this);\r
 224   }\r
 225   if(!fClusters) return;\r
 226   if(!fDetFuncList){\r
 227     AliDebug(1, "No detector functor list defined for task");\r
 228     return;\r
 229   }\r
 230   TIter detIter(fDetFuncList);\r
 231   for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){\r
 232     if(!(fkClusters = (TObjArray*)fClusters->At(idet))) continue;\r
 233     TMethodCall *det(NULL);\r
 234     detIter.Reset();\r
 235     while((det=dynamic_cast<TMethodCall*>(detIter()))) det->Execute(this);\r
 236   }\r
 237 }\r
 238 \r
 239 //_______________________________________________________\r
 240 Bool_t AliTRDrecoTask::GetRefFigure(Int_t /*ifig*/)\r
 241 {\r
 242   AliWarning("Retrieving reference figures not implemented.");\r
 243   return kFALSE;\r
 244 }\r
 245 \r
 246 //_______________________________________________________\r
 247 Bool_t AliTRDrecoTask::PutTrendValue(const Char_t *name, Double_t val, Double_t err)\r
 248 {\r
 249 // Generic publisher for trend values\r
 250 \r
 251   AliTRDtrendingManager *tm = AliTRDtrendingManager::Instance();\r
 252   tm->AddValue(Form("%s_%s", GetName(), name), val, err);\r
 253   return kTRUE;\r
 254 }\r
 255 \r
 256 //_______________________________________________________\r
 257 void AliTRDrecoTask::InitFunctorList()\r
 258 {\r
 259 // Initialize list of functors\r
 260 \r
 261   TClass *c = this->IsA();\r
 262   if(fPlotFuncList) fPlotFuncList->Clear();\r
 263   if(fDetFuncList) fDetFuncList->Clear();\r
 264 \r
 265   TMethod *m(NULL);\r
 266   TIter methIter(c->GetListOfMethods());\r
 267   while((m=dynamic_cast<TMethod*>(methIter()))){\r
 268     TString name(m->GetName());\r
 269     if(name.BeginsWith("Plot")){\r
 270       if(!fPlotFuncList) fPlotFuncList = new TList();\r
 271       fPlotFuncList->AddLast(new TMethodCall(c, (const char*)name, ""));\r
 272     } else if(name.BeginsWith("Det")){\r
 273       if(!fDetFuncList) fDetFuncList = new TList();\r
 274       fDetFuncList->AddLast(new TMethodCall(c, (const char*)name, ""));\r
 275     }\r
 276   }\r
 277 }\r
 278 \r
 279 //_______________________________________________________\r
 280 Bool_t AliTRDrecoTask::Load(const Char_t *file, const Char_t *dir)\r
 281 {\r
 282 // Generic container loader\r
 283 \r
 284   if(!TFile::Open(file)){\r
 285     AliWarning(Form("Couldn't open file %s.", file));\r
 286     return kFALSE;\r
 287   }\r
 288   if(!gFile->cd(dir)){\r
 289     AliWarning(Form("Couldn't cd to %s in %s.", dir, file));\r
 290     return kFALSE;\r
 291   }\r
 292   TObjArray *o = NULL;\r
 293   if(!(o = (TObjArray*)gDirectory->Get(GetName()))){\r
 294     AliWarning("Missing histogram container.");\r
 295     return kFALSE;\r
 296   }\r
 297   fContainer = (TObjArray*)o->Clone(GetName());\r
 298   gFile->Close();\r
 299   return kTRUE;\r
 300 }\r
 301 \r
 302 //________________________________________________________\r
 303 Bool_t AliTRDrecoTask::LoadDetectorMap(const Char_t *file, const Char_t *dir)\r
 304 {\r
 305 // Load detector map.\r
 306 \r
 307   if(!TFile::Open(file)){\r
 308     AliWarning(Form("Couldn't open file %s.", file));\r
 309     return kFALSE;\r
 310   }\r
 311   if(!gFile->cd(dir)){\r
 312     AliWarning(Form("Couldn't cd to %s in %s.", dir, file));\r
 313     return kFALSE;\r
 314   }\r
 315   TObjArray *info = NULL;\r
 316   if(!(info = (TObjArray*)gDirectory->Get("TRDinfoGen"))){\r
 317     AliWarning("Missing TRDinfoGen container.");\r
 318     return kFALSE;\r
 319   }\r
 320   TObjArray *dets = (TObjArray*)info->FindObject("Chambers");\r
 321   if(!dets){\r
 322     TVector *vdets = (TVector*)info->At(4);\r
 323     if(!vdets){\r
 324       AliWarning("Missing detector map from TRDinfoGen results.");\r
 325       return kFALSE;\r
 326     } else fDetsV = (TVector*)vdets->Clone();\r
 327   } else fDets = (TObjArray*)dets->Clone("Chambers");\r
 328   gFile->Close();\r
 329   return kTRUE;\r
 330 }\r
 331 \r
 332 \r
 333 //________________________________________________________\r
 334 Bool_t AliTRDrecoTask::Save(TObjArray * const results){\r
 335   //\r
 336   // Store the output graphs in a ROOT file\r
 337   // Input TObject array will not be written as Key to the file,\r
 338   // only content itself\r
 339   //\r
 340 \r
 341   TDirectory *cwd = gDirectory;\r
 342   if(!TFile::Open(Form("TRD.Result%s.root", GetName()), "RECREATE")) return kFALSE;\r
 343 \r
 344   TIterator *iter = results->MakeIterator();\r
 345   TObject *inObject = NULL, *outObject = NULL;\r
 346   while((inObject = iter->Next())){\r
 347     outObject = inObject->Clone();\r
 348     outObject->Write(NULL, TObject::kSingleKey);\r
 349   }\r
 350   delete iter;\r
 351   gFile->Close(); delete gFile;\r
 352   cwd->cd(); \r
 353   return kTRUE;\r
 354 }\r
 355 \r
 356 //_______________________________________________________\r
 357 Bool_t AliTRDrecoTask::PostProcess()\r
 358 {\r
 359 // To be implemented by particular tasks\r
 360 \r
 361   AliWarning("Post processing of reference histograms not implemented.");\r
 362   return kTRUE;\r
 363 }\r
 364 \r
 365 //_______________________________________________________\r
 366 void AliTRDrecoTask::MakeDetectorPlot(Int_t ly, const Option_t *opt)\r
 367 {\r
 368 // Draw chamber boundaries in eta/phi plots with misalignments\r
 369 // based on info collected by AliTRDinfoGen\r
 370 \r
 371   if(!fDets){\r
 372     AliWarning("OLD Detector map and status not available. Try NEW");\r
 373     MakeDetectorPlotNEW(ly, opt);\r
 374     return;\r
 375   }\r
 376   Float_t xmin(0.), xmax(0.);\r
 377   TBox *gdet = new TBox();\r
 378   gdet->SetLineColor(kBlack);gdet->SetFillColor(kBlack);\r
 379   Int_t style[] = {0, 3003};\r
 380   for(Int_t idet(0); idet<540; idet++){\r
 381     if(idet%6 != ly) continue;\r
 382     TVectorF *det((TVectorF*)fDets->At(idet));\r
 383     if(!det) continue;\r
 384     Int_t iopt = Int_t((*det)[4]);\r
 385     if(strcmp(opt, "eta")==0){\r
 386       xmin=(*det)[0]; xmax=(*det)[2];\r
 387     } else if(strcmp(opt, "pad")==0){\r
 388       Int_t stk(AliTRDgeometry::GetStack(idet));\r
 389       xmin=-0.6+16*(4-stk)-(stk<2?4:0); xmax=xmin+(stk==2?12:16)-0.2;\r
 390     } else continue;\r
 391     AliDebug(2, Form("det[%03d] 0[%+4.1f(%+4.1f) %+4.1f] 1[%+4.1f(%+4.1f) %+4.1f] opt[%d]", idet, xmin, (*det)[0], (*det)[1], xmax, (*det)[2], (*det)[3], iopt));\r
 392     if(iopt==1){\r
 393       gdet->SetFillStyle(style[1]);gdet->SetFillColor(kBlack);\r
 394       gdet->DrawBox(xmin, (*det)[1], xmax, (*det)[3]);\r
 395     } else {\r
 396       gdet->SetFillStyle(style[0]);\r
 397       gdet->DrawBox(xmin, (*det)[1], xmax, (*det)[3]);\r
 398       if(iopt==2){\r
 399         gdet->SetFillStyle(style[1]);gdet->SetFillColor(kGreen);\r
 400         gdet->DrawBox(xmin, (*det)[1], xmax, 0.5*((*det)[3]+(*det)[1]));\r
 401       } else if(iopt==3){\r
 402         gdet->SetFillStyle(style[1]);gdet->SetFillColor(kRed);\r
 403         gdet->DrawBox(xmin, 0.5*((*det)[3]+(*det)[1]), xmax, (*det)[3]);\r
 404       } else if(iopt!=0) AliError(Form("Wrong chmb. status[%d] for det[%03d]", iopt, idet));\r
 405     }\r
 406   }\r
 407   Float_t dsm = TMath::TwoPi()/AliTRDgeometry::kNsector;\r
 408   xmin=0.;\r
 409   if(strcmp(opt, "pad")==0) xmin=38.;\r
 410   TLatex *sm = new TLatex(); sm->SetTextAlign(22);sm->SetTextColor(0); sm->SetTextFont(32);sm->SetTextSize(0.03);\r
 411   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
 412 }\r
 413 \r
 414 //_______________________________________________________\r
 415 void AliTRDrecoTask::MakeDetectorPlotNEW(Int_t ly, const Option_t *opt)\r
 416 {\r
 417 // Draw chamber boundaries in eta/phi plots with misalignments\r
 418 // based on info collected by AliTRDinfoGen NEW data storage\r
 419 \r
 420   if(!fDetsV){\r
 421     AliWarning("NEW Detector map and status not available.");\r
 422     return;\r
 423   }\r
 424   Float_t xmin(0.), xmax(0.);\r
 425   TBox *gdet = new TBox();\r
 426   gdet->SetLineColor(kBlack);gdet->SetFillColor(kBlack);\r
 427   Int_t style[] = {0, 3003};\r
 428   for(Int_t idet(0), jdet(0); idet<AliTRDgeometry::kNdet; idet++, jdet+=5){\r
 429     if(idet%6 != ly) continue;\r
 430     Int_t iopt = Int_t((*fDetsV)[jdet+4]);\r
 431     if(strcmp(opt, "eta")==0){\r
 432       xmin=(*fDetsV)[jdet+0]; xmax=(*fDetsV)[jdet+2];\r
 433     } else if(strcmp(opt, "pad")==0){\r
 434       Int_t stk(AliTRDgeometry::GetStack(idet));\r
 435       xmin=-0.6+16*(4-stk)-(stk<2?4:0); xmax=xmin+(stk==2?12:16)-0.2;\r
 436     } else continue;\r
 437     AliDebug(2, Form("det[%03d] 0[%+4.1f(%+4.1f) %+4.1f] 1[%+4.1f(%+4.1f) %+4.1f] opt[%d]", idet, xmin, (*fDetsV)[jdet+0], (*fDetsV)[jdet+1], xmax, (*fDetsV)[jdet+2], (*fDetsV)[jdet+3], iopt));\r
 438     if(iopt==1){\r
 439       gdet->SetFillStyle(style[1]);gdet->SetFillColor(kBlack);\r
 440       gdet->DrawBox(xmin, (*fDetsV)[jdet+1], xmax, (*fDetsV)[jdet+3]);\r
 441     } else {\r
 442       gdet->SetFillStyle(style[0]);\r
 443       gdet->DrawBox(xmin, (*fDetsV)[jdet+1], xmax, (*fDetsV)[jdet+3]);\r
 444       if(iopt==2){\r
 445         gdet->SetFillStyle(style[1]);gdet->SetFillColor(kGreen);\r
 446         gdet->DrawBox(xmin, (*fDetsV)[jdet+1], xmax, 0.5*((*fDetsV)[jdet+3]+(*fDetsV)[jdet+1]));\r
 447       } else if(iopt==3){\r
 448         gdet->SetFillStyle(style[1]);gdet->SetFillColor(kRed);\r
 449         gdet->DrawBox(xmin, 0.5*((*fDetsV)[jdet+3]+(*fDetsV)[jdet+1]), xmax, (*fDetsV)[jdet+3]);\r
 450       } else if(iopt!=0) AliError(Form("Wrong chmb. status[%d] for det[%03d]", iopt, idet));\r
 451     }\r
 452   }\r
 453   Float_t dsm = TMath::TwoPi()/AliTRDgeometry::kNsector;\r
 454   xmin=0.;\r
 455   if(strcmp(opt, "pad")==0) xmin=38.;\r
 456   TLatex *sm = new TLatex(); sm->SetTextAlign(22);sm->SetTextColor(0); sm->SetTextFont(32);sm->SetTextSize(0.03);\r
 457   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
 458 }\r
 459 \r
 460 \r
 461 //_______________________________________________________\r
 462 void AliTRDrecoTask::MakeSummary()\r
 463 {\r
 464 // To be implemented by particular tasks\r
 465   AliWarning("Summary not available");\r
 466 }\r
 467 \r
 468 //_______________________________________________________\r
 469 void AliTRDrecoTask::SetDebugLevel(Int_t level)\r
 470 {\r
 471 // Generic debug handler\r
 472 \r
 473   AliAnalysisTaskSE::SetDebugLevel(level);\r
 474   if(DebugLevel()>=1){\r
 475     AliInfo(Form("Debug Level for Task %s set to %d", GetName(), level));\r
 476     TDirectory *savedir = gDirectory;\r
 477     fgDebugStream = new TTreeSRedirector("TRD.DebugPerformance.root");\r
 478     savedir->cd();\r
 479   }\r
 480 }\r
 481 \r
 482 //____________________________________________________________________\r
 483 void AliTRDrecoTask::Terminate(Option_t *)\r
 484 {\r
 485   //\r
 486   // Terminate\r
 487   //\r
 488 \r
 489   if(fgDebugStream){ \r
 490     delete fgDebugStream;\r
 491     fgDebugStream = NULL;\r
 492   }\r
 493   fContainer = dynamic_cast<TObjArray *>(GetOutputData(1));\r
 494   if(fContainer && fRunTerminate){\r
 495     PostProcess();\r
 496     MakeSummary();\r
 497   }\r
 498 }\r
 499 \r
 500 //________________________________________________________\r
 501 void AliTRDrecoTask::SetNormZ(TH2 *h2, Int_t bxmin, Int_t bxmax, Int_t bymin, Int_t bymax, Float_t thr)\r
 502 {\r
 503 // Normalize histo content to the mean value in the range specified by bin ranges\r
 504 // [bxmin, bxmax] on the x axis and [bymin, bymax] on the y axis.\r
 505 // Optionally a threshold "thr" can be specified to disregard entries with no meaning\r
 506 \r
 507   Float_t s = 0., c=0.; Int_t is(0);\r
 508   for(Int_t ix(bxmin); ix<=(bxmax>0?bxmax:(h2->GetXaxis()->GetNbins())); ix++){\r
 509     for(Int_t iy(bymin); iy<=(bymax>0?bymax:(h2->GetYaxis()->GetNbins())); iy++){\r
 510       if((c = h2->GetBinContent(ix, iy))<thr) continue;\r
 511       s += c; is++;\r
 512     }\r
 513   }\r
 514   s/= (is?is:1);\r
 515   for(Int_t ix(1); ix<=h2->GetXaxis()->GetNbins(); ix++){\r
 516     for(Int_t iy(1); iy<=h2->GetYaxis()->GetNbins(); iy++){\r
 517       if((c = h2->GetBinContent(ix, iy))<thr) h2->SetBinContent(ix, iy, thr-1000);\r
 518       else h2->SetBinContent(ix, iy, 100.*(c/s-1.));\r
 519     }\r
 520   }\r
 521 }\r
 522 \r
 523 //________________________________________________________\r
 524 void AliTRDrecoTask::SetRangeZ(TH2 *h2, Float_t min, Float_t max, Float_t thr)\r
 525 {\r
 526 // Set range on Z axis such to avoid outliers\r
 527 \r
 528   Float_t c(0.), dz(1.e-3*(max-min));\r
 529   for(Int_t ix(1); ix<=h2->GetXaxis()->GetNbins(); ix++){\r
 530     for(Int_t iy(1); iy<=h2->GetYaxis()->GetNbins(); iy++){\r
 531       if((c = h2->GetBinContent(ix, iy))<thr) continue;\r
 532       if(c<=min) h2->SetBinContent(ix, iy, min+dz);\r
 533     }\r
 534   }\r
 535   h2->GetZaxis()->SetRangeUser(min, max);\r
 536 }\r
 537 \r
 538 //________________________________________________________\r
 539 Float_t AliTRDrecoTask::GetMeanStat(TH1 *h, Float_t cut, Option_t *opt)\r
 540 {\r
 541 // return mean number of entries/bin of histogram "h"\r
 542 // if option "opt" is given the following values are accepted:\r
 543 //   "<" : consider only entries less than "cut"\r
 544 //   ">" : consider only entries greater than "cut"\r
 545 \r
 546   //Int_t dim(h->GetDimension());\r
 547   Int_t nbx(h->GetNbinsX()), nby(h->GetNbinsY()), nbz(h->GetNbinsZ());\r
 548   Double_t sum(0.); Int_t n(0);\r
 549   for(Int_t ix(1); ix<=nbx; ix++)\r
 550     for(Int_t iy(1); iy<=nby; iy++)\r
 551       for(Int_t iz(1); iz<=nbz; iz++){\r
 552         if(strcmp(opt, "")==0){sum += h->GetBinContent(ix, iy, iz); n++;}\r
 553         else{\r
 554           if(strcmp(opt, "<")==0) {\r
 555             if(h->GetBinContent(ix, iy, iz)<cut) {sum += h->GetBinContent(ix, iy, iz); n++;}\r
 556           } else if(strcmp(opt, ">")==0){\r
 557             if(h->GetBinContent(ix, iy, iz)>cut) {sum += h->GetBinContent(ix, iy, iz); n++;}\r
 558           } else {sum += h->GetBinContent(ix, iy, iz); n++;}\r
 559         }\r
 560       }\r
 561   return n>0?sum/n:0.;\r
 562 }\r
 563 \r
 564 //________________________________________________________\r
 565 AliTRDrecoTask::AliTRDrecoProjection::AliTRDrecoProjection()\r
 566   :TNamed()\r
 567   ,fH(NULL)\r
 568   ,fNrebin(0)\r
 569   ,fRebinX(NULL)\r
 570   ,fRebinY(NULL)\r
 571 {\r
 572   // constructor\r
 573   memset(fAx, 0, 3*sizeof(Int_t));\r
 574   memset(fRange, 0, 4*sizeof(Float_t));\r
 575 }\r
 576 \r
 577 //________________________________________________________\r
 578 AliTRDrecoTask::AliTRDrecoProjection::~AliTRDrecoProjection()\r
 579 {\r
 580   // destructor\r
 581   if(fH) delete fH;\r
 582 }\r
 583 \r
 584 //________________________________________________________\r
 585 void AliTRDrecoTask::AliTRDrecoProjection::Build(const Char_t *n, const Char_t *t, Int_t ix, Int_t iy, Int_t iz, TAxis *aa[])\r
 586 {\r
 587 // check and build (if neccessary) projection determined by axis "ix", "iy" and "iz"\r
 588   if(!aa[ix] || !aa[iy] || !aa[iz]) return;\r
 589   TAxis *ax(aa[ix]), *ay(aa[iy]), *az(aa[iz]);\r
 590   // check ax definiton to protect against older versions of the data\r
 591   if(ax->GetNbins()<=0 || (ax->GetXmax()-ax->GetXmin())<=0.){\r
 592     AliWarning(Form("Wrong definition of axis[%d] \"%s\"[%d](%f %f).", ix, ax->GetTitle(), ax->GetNbins(), ax->GetXmin(), ax->GetXmax()));\r
 593     return;\r
 594   }\r
 595   if(ay->GetNbins()<=0 || (ay->GetXmax()-ay->GetXmin())<=0.){\r
 596     AliWarning(Form("Wrong definition of axis[%d] \"%s\"[%d](%f %f).", ix, ay->GetTitle(), ay->GetNbins(), ay->GetXmin(), ay->GetXmax()));\r
 597     return;\r
 598   }\r
 599   if(az->GetNbins()<=0 || (az->GetXmax()-az->GetXmin())<=0.){\r
 600     AliWarning(Form("Wrong definition of axis[%d] \"%s\"[%d](%f %f).", ix, az->GetTitle(), az->GetNbins(), az->GetXmin(), az->GetXmax()));\r
 601     return;\r
 602   }\r
 603   SetNameTitle(n,t);\r
 604   fH = new TH3I(n, Form("%s;%s;%s;%s", t, ax->GetTitle(), ay->GetTitle(), az->GetTitle()),\r
 605     ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),\r
 606     ay->GetNbins(), ay->GetXmin(), ay->GetXmax(),\r
 607     az->GetNbins(), az->GetXmin(), az->GetXmax());\r
 608   fAx[0] = ix; fAx[1] = iy; fAx[2] = iz;\r
 609   fRange[0] = az->GetXmin()/3.; fRange[1] = az->GetXmax()/3.;\r
 610   AliDebug(2, Form("H3(%s, %s) :: %s[%3d %4.2f %4.2f]%s[%3d %4.2f %4.2f]%s[%3d %4.2f %4.2f]", n, t,\r
 611     ax->GetTitle(), ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),\r
 612     ay->GetTitle(), ay->GetNbins(), ay->GetXmin(), ay->GetXmax(),\r
 613     az->GetTitle(), az->GetNbins(), az->GetXmin(), az->GetXmax()));\r
 614 }\r
 615 \r
 616 //________________________________________________________\r
 617 AliTRDrecoTask::AliTRDrecoProjection& AliTRDrecoTask::AliTRDrecoProjection::operator=(const AliTRDrecoProjection& rhs)\r
 618 {\r
 619 // copy projections\r
 620   if(this == &rhs) return *this;\r
 621 \r
 622   TNamed::operator=(rhs);\r
 623   if(fNrebin){fNrebin=0; delete [] fRebinX; delete [] fRebinY;}\r
 624   if(rhs.fNrebin) SetRebinStrategy(rhs.fNrebin, rhs.fRebinX, rhs.fRebinY);\r
 625   memcpy(fAx, rhs.fAx, 3*sizeof(Int_t));\r
 626   memcpy(fRange, rhs.fRange, 4*sizeof(Float_t));\r
 627   if(fH) delete fH;\r
 628   if(rhs.fH) fH=(TH3I*)rhs.fH->Clone(Form("%s_CLONE", rhs.fH->GetName()));\r
 629   return *this;\r
 630 }\r
 631 \r
 632 //________________________________________________________\r
 633 AliTRDrecoTask::AliTRDrecoProjection& AliTRDrecoTask::AliTRDrecoProjection::operator+=(const AliTRDrecoProjection& other)\r
 634 {\r
 635 // increment projections\r
 636   if(!fH || !other.fH) return *this;\r
 637   AliDebug(2, Form("%s+=%s [%s+=%s]", GetName(), other.GetName(), fH->GetName(), (other.fH)->GetName()));\r
 638   fH->Add(other.fH);\r
 639   return *this;\r
 640 }\r
 641 \r
 642 //________________________________________________________\r
 643 void AliTRDrecoTask::AliTRDrecoProjection::Increment(Int_t bin[], Double_t v)\r
 644 {\r
 645 // increment bin with value "v" pointed by general coord in "bin"\r
 646   if(!fH) return;\r
 647   AliDebug(4, Form("  %s[%2d]", fH->GetName(), Int_t(v)));\r
 648   fH->AddBinContent(fH->GetBin(bin[fAx[0]],bin[fAx[1]],bin[fAx[2]]), Int_t(v));\r
 649 }\r
 650 \r
 651 //________________________________________________________\r
 652 TH2* AliTRDrecoTask::AliTRDrecoProjection::Projection2D(const Int_t nstat, const Int_t ncol, const Int_t mid, Bool_t del)\r
 653 {\r
 654 // build the 2D projection and adjust binning\r
 655 \r
 656   const Char_t *title[] = {"Mean", "#mu", "MPV"};\r
 657   if(!fH){\r
 658     AliDebug(1, Form("Missing 3D in %s", GetName()));\r
 659     return NULL;\r
 660   }\r
 661   TAxis *ax(fH->GetXaxis()), *ay(fH->GetYaxis()), *az(fH->GetZaxis());\r
 662   TH2D *h2s(NULL), *hyx(NULL);\r
 663   if(!(h2s = (TH2D*)fH->Project3D("yx"))){\r
 664     AliDebug(1, Form("Failed Project3D(\"yx\") in %s", GetName()));\r
 665     return NULL;\r
 666   }\r
 667   // save a copy of the original distribution\r
 668   if(!del){\r
 669     hyx = (TH2D*)h2s->Clone();\r
 670     hyx->SetName(Form("%sEn", fH->GetName()));\r
 671   }\r
 672   Int_t irebin(0), dxBin(1), dyBin(1);\r
 673   while(irebin<fNrebin && (AliTRDrecoTask::GetMeanStat(h2s, .5, ">")<nstat)){\r
 674     h2s->Rebin2D(fRebinX[irebin], fRebinY[irebin]);\r
 675     dxBin*=fRebinX[irebin];dyBin*=fRebinY[irebin];\r
 676     irebin++;\r
 677   }\r
 678   Int_t nx(h2s->GetNbinsX()), ny(h2s->GetNbinsY());\r
 679   delete h2s;\r
 680   if(mid<0) return NULL;\r
 681 \r
 682   // start projection\r
 683   TH1 *h(NULL); Int_t n(0);\r
 684   Float_t dz=(fRange[1]-fRange[1])/ncol;\r
 685   TString titlez(az->GetTitle()); TObjArray *tokenTitle(titlez.Tokenize(" "));\r
 686   Int_t nt(tokenTitle->GetEntriesFast());\r
 687   TH2 *h2(NULL);\r
 688   if((h2 = (TH2*)gDirectory->Get(Form("%s_2D", fH->GetName())))) delete h2; // avoid ROOT warning messages\r
 689   h2 = new TH2F(Form("%s_2D", fH->GetName()),\r
 690             Form("%s;%s;%s;%s(%s) %s", fH->GetTitle(), ax->GetTitle(), ay->GetTitle(), title[mid], nt>0?(*tokenTitle)[0]->GetName():"", nt>1?(*tokenTitle)[1]->GetName():""),\r
 691             nx, ax->GetXmin(), ax->GetXmax(), ny, ay->GetXmin(), ay->GetXmax());\r
 692   h2->SetContour(ncol);\r
 693   h2->GetZaxis()->CenterTitle();\r
 694   h2->GetZaxis()->SetTitleOffset(1.4);\r
 695   h2->GetZaxis()->SetRangeUser(fRange[0], fRange[1]);\r
 696   AliDebug(2, Form("%s[%s] nx[%d] ny[%d]", h2->GetName(), h2->GetTitle(), nx, ny));\r
 697   for(Int_t iy(0); iy<ny; iy++){\r
 698     for(Int_t ix(0); ix<nx; ix++){\r
 699       h = fH->ProjectionZ(Form("%s_z", h2->GetName()), ix*dxBin+1, (ix+1)*dxBin, iy*dyBin+1, (iy+1)*dyBin);\r
 700       Int_t ne((Int_t)h->Integral());\r
 701       //printf("  x[%2d %2d] y[%2d %2d] ne[%4d]\n", ix*dxBin+1, (ix+1)*dxBin, iy*dyBin+1, (iy+1)*dyBin, ne);\r
 702       if(ne<nstat/2){\r
 703         h2->SetBinContent(ix+1, iy+1, -999);\r
 704         h2->SetBinError(ix+1, iy+1, 1.);\r
 705         n++;\r
 706       }else{\r
 707         // redo the projection by adding 1 bin @ left and 1 bin @ right for smoothing\r
 708         h = fH->ProjectionZ(Form("%s_z", h2->GetName()), ix*dxBin, (ix+1)*dxBin+1, iy*dyBin, (iy+1)*dyBin+1);\r
 709         Float_t v(h->GetMean()), ve(h->GetRMS());\r
 710         if(mid==1){\r
 711           TF1 fg("fg", "gaus", az->GetXmin(), az->GetXmax());\r
 712           fg.SetParameter(0, Float_t(ne)); fg.SetParameter(1, v); fg.SetParameter(2, ve);\r
 713           h->Fit(&fg, "WQ0");\r
 714           v = fg.GetParameter(1); ve = fg.GetParameter(2);\r
 715         } else if (mid==2) {\r
 716           TF1 fl("fl", "landau", az->GetXmin(), az->GetXmax());\r
 717           fl.SetParameter(0, Float_t(ne)); fl.SetParameter(1, v); fl.SetParameter(2, ve);\r
 718           h->Fit(&fl, "WQ0");\r
 719           v = fl.GetMaximumX(); ve = fl.GetParameter(2);\r
 720 /*          TF1 fgle("gle", "[0]*TMath::Landau(x, [1], [2], 1)*TMath::Exp(-[3]*x/[1])", az->GetXmin(), az->GetXmax());\r
 721           fgle.SetParameter(0, fl.GetParameter(0));\r
 722           fgle.SetParameter(1, fl.GetParameter(1));\r
 723           fgle.SetParameter(2, fl.GetParameter(2));\r
 724           fgle.SetParameter(3, 1.);fgle.SetParLimits(3, 0., 5.);\r
 725           h->Fit(&fgle, "WQ");\r
 726           v = fgle.GetMaximumX(); ve = fgle.GetParameter(2);*/\r
 727         }\r
 728         if(v<fRange[0]) h2->SetBinContent(ix+1, iy+1, fRange[0]+0.1*dz);\r
 729         else h2->SetBinContent(ix+1, iy+1, v);\r
 730         h2->SetBinError(ix+1, iy+1, ve);\r
 731       }\r
 732     }\r
 733   }\r
 734   if(h) delete h;\r
 735   if(n==nx*ny){  // clean empty projections\r
 736     AliDebug(1, Form("Empty projection in %s", GetName()));\r
 737     delete h2; h2=NULL;\r
 738   }\r
 739   return h2;\r
 740 }\r
 741 \r
 742 //________________________________________________________\r
 743 void AliTRDrecoTask::AliTRDrecoProjection::SetRebinStrategy(Int_t n, Int_t rebx[], Int_t reby[])\r
 744 {\r
 745 // define rebinning strategy for this projection\r
 746   fNrebin = n;\r
 747   fRebinX = new Int_t[n]; memcpy(fRebinX, rebx, n*sizeof(Int_t));\r
 748   fRebinY = new Int_t[n]; memcpy(fRebinY, reby, n*sizeof(Int_t));\r
 749 }\r
 750 \r
 751 \r
 752 \r