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