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/AliTRDchmbInfo.h"\r
  41 #include "info/AliTRDeventInfo.h"\r
  42 #include "info/AliTRDtrendingManager.h"\r
  43 #include "AliTRDrecoTask.h"\r
  44 #include "AliTRDtrackV1.h"\r
  45 #include "AliTRDpidUtil.h"\r
  46 \r
  47 ClassImp(AliTRDrecoTask)\r
  48 \r
  49 Float_t AliTRDrecoTask::fgPt[AliTRDrecoTask::fgNPt+1] = {0.};\r
  50 TTreeSRedirector* AliTRDrecoTask::fgDebugStream(NULL);\r
  51 TH1* AliTRDrecoTask::fgProjector(NULL);\r
  52 //_______________________________________________________\r
  53 AliTRDrecoTask::AliTRDrecoTask()\r
  54   : AliAnalysisTaskSE()\r
  55   ,fNRefFigures(0)\r
  56   ,fDets(NULL)\r
  57   ,fDetsV(NULL)\r
  58   ,fContainer(NULL)\r
  59   ,fEvent(NULL)\r
  60   ,fTracks(NULL)\r
  61   ,fClusters(NULL)\r
  62   ,fkClusters(NULL)\r
  63   ,fkTrack(NULL)\r
  64   ,fkMC(NULL)\r
  65   ,fkESD(NULL)\r
  66   ,fSpecies(-6)\r
  67   ,fTriggerSlot(0)\r
  68   ,fPt(-1.)\r
  69   ,fPhi(0.)\r
  70   ,fEta(0.)\r
  71   ,fNpt(0)\r
  72   ,fTriggerList(NULL)\r
  73   ,fPlotFuncList(NULL)\r
  74   ,fDetFuncList(NULL)\r
  75   ,fRunTerminate(kFALSE)\r
  76 {\r
  77 // Default constructor\r
  78   snprintf(fNameId, 10, "no name");\r
  79 }\r
  80 \r
  81 //_______________________________________________________\r
  82 AliTRDrecoTask::AliTRDrecoTask(const char *name, const char *title)\r
  83   : AliAnalysisTaskSE(name)\r
  84   ,fNRefFigures(0)\r
  85   ,fDets(NULL)\r
  86   ,fDetsV(NULL)\r
  87   ,fContainer(NULL)\r
  88   ,fEvent(NULL)\r
  89   ,fTracks(NULL)\r
  90   ,fClusters(NULL)\r
  91   ,fkClusters(NULL)\r
  92   ,fkTrack(NULL)\r
  93   ,fkMC(NULL)\r
  94   ,fkESD(NULL)\r
  95   ,fSpecies(-6)\r
  96   ,fTriggerSlot(0)\r
  97   ,fPt(-1.)\r
  98   ,fPhi(0.)\r
  99   ,fEta(0.)\r
 100   ,fNpt(0)\r
 101   ,fTriggerList(NULL)\r
 102   ,fPlotFuncList(NULL)\r
 103   ,fDetFuncList(NULL)\r
 104   ,fRunTerminate(kFALSE)\r
 105 {\r
 106 // Constructor for all derived performance tasks\r
 107 \r
 108   SetTitle(title);\r
 109   snprintf(fNameId, 10, "no name");\r
 110   DefineInput (1, TObjArray::Class()); // track list\r
 111   DefineInput (2, AliTRDeventInfo::Class()); // event info object\r
 112   DefineInput (3, TObjArray::Class()); // cluster list object\r
 113   DefineOutput(1, TObjArray::Class()); // histogram list\r
 114 }\r
 115 \r
 116 //_______________________________________________________\r
 117 AliTRDrecoTask::~AliTRDrecoTask() \r
 118 {\r
 119 \r
 120   // Generic task destructor\r
 121 \r
 122   AliDebug(2, Form(" Ending task %s[%s]", GetName(), GetTitle()));\r
 123  if(fgDebugStream){\r
 124     delete fgDebugStream;\r
 125     fgDebugStream = NULL;\r
 126   }\r
 127 \r
 128   if(fgProjector){\r
 129     delete fgProjector;\r
 130     fgProjector = NULL;\r
 131   }\r
 132 \r
 133   if(fPlotFuncList){\r
 134     fPlotFuncList->Delete();\r
 135     delete fPlotFuncList;\r
 136     fPlotFuncList = NULL;\r
 137   }\r
 138   if(fDetFuncList){\r
 139     fDetFuncList->Delete();\r
 140     delete fDetFuncList;\r
 141     fDetFuncList = NULL;\r
 142   }\r
 143   \r
 144   if(fDets){\r
 145     if(fDets->IsOwner()) fDets->Delete();\r
 146     delete fDets;\r
 147     fDets = NULL;\r
 148   }\r
 149   if(fDetsV) delete fDetsV; fDetsV=NULL;\r
 150   if(fTriggerList){fTriggerList->Delete(); delete fTriggerList;}\r
 151 \r
 152   if(fContainer && !(AliAnalysisManager::GetAnalysisManager() && AliAnalysisManager::GetAnalysisManager()->IsProofMode())){\r
 153     if(fContainer->IsOwner()) fContainer->Delete();\r
 154     delete fContainer;\r
 155     fContainer = NULL;\r
 156   }\r
 157 \r
 158 /*  if(fgTrendPoint){\r
 159     TFile::Open("TRD.PerformanceTrend.root", "UPDATE");\r
 160     fgTrendPoint->Write();\r
 161     delete fgTrendPoint;\r
 162     fgTrendPoint=NULL;\r
 163     gFile->Close();\r
 164   }*/\r
 165 }\r
 166 \r
 167 //_______________________________________________________\r
 168 Int_t AliTRDrecoTask::GetNRefFigures() const  \r
 169 { \r
 170   if(!fNRefFigures) AliWarning("No reference plots available.");\r
 171   return fNRefFigures; \r
 172 } \r
 173 \r
 174 //____________________________________________________________________\r
 175 Int_t AliTRDrecoTask::GetPtBin(Float_t pt)\r
 176 {\r
 177 // Get significant (very low, low, medium, high, very high) pt bin\r
 178 \r
 179   Int_t ipt(0);\r
 180   while(ipt<fNpt){\r
 181     if(pt<fgPt[ipt]) break;\r
 182     ipt++;\r
 183   }\r
 184   return ipt-1;\r
 185 }\r
 186 \r
 187 //_______________________________________________________\r
 188 Bool_t AliTRDrecoTask::MakeMomSegmentation()\r
 189 {\r
 190   switch(fNpt){\r
 191   case fgNPt:\r
 192     fgPt[0]=0.3;\r
 193     for(Int_t j(1); j<=fgNPt; j++) fgPt[j]=fgPt[j-1]+(TMath::Exp(j*j*2.e-3)-1.);\r
 194     AliDebug(2, "Using debug momentum segmentation");\r
 195     break;\r
 196   case 4:\r
 197     fgPt[0]=0.5; fgPt[1]=0.8; fgPt[2]=1.5; fgPt[3]=5.;\r
 198     AliDebug(2, "Using default momentum segmentation");\r
 199     break;\r
 200   default:\r
 201     AliError(Form("Momentum segmentation %d not supported.", fNpt));\r
 202     fNpt=0;\r
 203     return kFALSE;\r
 204   }\r
 205   return kTRUE;\r
 206 }\r
 207 \r
 208 //_______________________________________________________\r
 209 void AliTRDrecoTask::UserCreateOutputObjects()\r
 210 {\r
 211   if(!HasFunctorList()) InitFunctorList();\r
 212   if(DebugLevel()) fNpt = fgNPt;\r
 213   else fNpt = 4;\r
 214   MakeMomSegmentation();\r
 215 \r
 216   fContainer = Histos();\r
 217   PostData(1, fContainer);\r
 218 }\r
 219 \r
 220 //_______________________________________________________\r
 221 void AliTRDrecoTask::UserExec(Option_t *)\r
 222 {\r
 223 // Loop over Plot functors published by particular tasks\r
 224 \r
 225   fTracks   = dynamic_cast<TObjArray *>(GetInputData(1));\r
 226   fEvent    = dynamic_cast<AliTRDeventInfo *>(GetInputData(2));\r
 227   fTriggerSlot=0;\r
 228   if(fTriggerList){\r
 229     for(Int_t itrig(0); itrig<fTriggerList->GetEntries(); itrig++){\r
 230       if(!fEvent->GetFiredTriggerClasses().Contains(((TObjString*)(*fTriggerList)[itrig])->GetName())) continue;\r
 231       //printf("\"%s\" selected\n", ((TObjString*)(*fTriggerList)[itrig])->GetName());\r
 232       SETBIT(fTriggerSlot,itrig);\r
 233     }\r
 234     if(!fTriggerSlot){\r
 235       AliDebug(2, Form("Triggers[%s] not used for %s", fEvent->GetFiredTriggerClasses().Data(),  GetName()));\r
 236       return;\r
 237     }\r
 238   }\r
 239   fClusters = dynamic_cast<TObjArray*>(GetInputData(3));\r
 240 \r
 241   if(!fPlotFuncList){\r
 242     AliWarning("No track functor list defined for the task");\r
 243     return;\r
 244   }\r
 245   if(!fTracks) return;\r
 246   if(!fTracks->GetEntriesFast()) return;\r
 247   else AliDebug(2, Form("Tracks[%d] for %s", fTracks->GetEntriesFast(), GetName()));\r
 248 \r
 249   Int_t itrk(-1);\r
 250   AliTRDtrackInfo *trackInfo(NULL);\r
 251   TIter plotIter(fPlotFuncList);\r
 252   TObjArrayIter trackIter(fTracks);\r
 253   while((trackInfo = dynamic_cast<AliTRDtrackInfo*>(trackIter()))){\r
 254     itrk++; fPt=-1; fEta=0.; fPhi=0.; fSpecies=-6;\r
 255     fkMC    = trackInfo->GetMCinfo();\r
 256     fkESD   = trackInfo->GetESDinfo();\r
 257     if((fkTrack = trackInfo->GetTrack())){\r
 258       // cache properties of the track at TRD entrance\r
 259       // check input track status\r
 260       AliExternalTrackParam *tin(NULL);\r
 261       if(!(tin = fkTrack->GetTrackIn())) AliDebug(2, Form("Missing TRD track[%d] :: entry point.", itrk));\r
 262       else {\r
 263         fPt   = tin->Pt();\r
 264         fEta  = tin->Eta();\r
 265         Double_t xyz[3];\r
 266         if(!tin->GetXYZ(xyz)) AliDebug(2, Form("Failed TRD track[%d] :: global track postion", itrk));\r
 267         else fPhi  = TMath::ATan2(xyz[1], xyz[0]);\r
 268         fSpecies= fkTrack->Charge()*(AliTRDpidUtil::Mass2Pid(fkTrack->GetMass())+1);\r
 269       }\r
 270     } else AliDebug(2, Form("Missing TRD track[%d].", itrk));\r
 271 \r
 272     TMethodCall *plot(NULL);\r
 273     plotIter.Reset();\r
 274     while((plot=dynamic_cast<TMethodCall*>(plotIter()))) plot->Execute(this);\r
 275   }\r
 276   if(!fClusters) return;\r
 277   if(!fDetFuncList){\r
 278     AliDebug(1, "No detector functor list defined for task");\r
 279     return;\r
 280   }\r
 281   TIter detIter(fDetFuncList);\r
 282   for(Int_t idet(0); idet<AliTRDgeometry::kNdet; idet++){\r
 283     if(!(fkClusters = (TObjArray*)fClusters->At(idet))) continue;\r
 284     TMethodCall *det(NULL);\r
 285     detIter.Reset();\r
 286     while((det=dynamic_cast<TMethodCall*>(detIter()))) det->Execute(this);\r
 287   }\r
 288 }\r
 289 \r
 290 //_______________________________________________________\r
 291 Bool_t AliTRDrecoTask::GetRefFigure(Int_t /*ifig*/)\r
 292 {\r
 293   AliWarning("Retrieving reference figures not implemented.");\r
 294   return kFALSE;\r
 295 }\r
 296 \r
 297 //_______________________________________________________\r
 298 Bool_t AliTRDrecoTask::PutTrendValue(const Char_t *name, Double_t val, Double_t err)\r
 299 {\r
 300 // Generic publisher for trend values\r
 301 \r
 302   AliTRDtrendingManager *tm = AliTRDtrendingManager::Instance();\r
 303   if(!tm){\r
 304     AliError("Wrong usage of the trending functionality. Could not instantiate AliTRDtrendingManager singleton.");\r
 305     return kFALSE;\r
 306   }\r
 307   tm->AddValue(Form("%s_%s", GetName(), name), val, err);\r
 308   return kTRUE;\r
 309 }\r
 310 \r
 311 //_______________________________________________________\r
 312 void AliTRDrecoTask::InitFunctorList()\r
 313 {\r
 314 // Initialize list of functors\r
 315 \r
 316   TClass *c = this->IsA();\r
 317   if(fPlotFuncList) fPlotFuncList->Clear();\r
 318   if(fDetFuncList) fDetFuncList->Clear();\r
 319 \r
 320   TMethod *m(NULL);\r
 321   TIter methIter(c->GetListOfMethods());\r
 322   while((m=dynamic_cast<TMethod*>(methIter()))){\r
 323     TString name(m->GetName());\r
 324     if(name.BeginsWith("Plot")){\r
 325       if(!fPlotFuncList) fPlotFuncList = new TList();\r
 326       fPlotFuncList->AddLast(new TMethodCall(c, (const char*)name, ""));\r
 327     } else if(name.BeginsWith("Det")){\r
 328       if(!fDetFuncList) fDetFuncList = new TList();\r
 329       fDetFuncList->AddLast(new TMethodCall(c, (const char*)name, ""));\r
 330     }\r
 331   }\r
 332 }\r
 333 \r
 334 //_______________________________________________________\r
 335 Bool_t AliTRDrecoTask::Load(const Char_t *file, const Char_t *dir)\r
 336 {\r
 337 // Generic container loader\r
 338 \r
 339   if(!TFile::Open(file)){\r
 340     AliWarning(Form("Couldn't open file %s.", file));\r
 341     return kFALSE;\r
 342   }\r
 343   if(!gFile->cd(dir)){\r
 344     AliWarning(Form("Couldn't cd to %s in %s.", dir, file));\r
 345     gFile->Close();\r
 346     return kFALSE;\r
 347   }\r
 348   if(!(fContainer = (TObjArray*)gDirectory->Get(GetName()))){\r
 349     AliWarning("Missing histogram container.");\r
 350     gFile->Close();\r
 351     return kFALSE;\r
 352   }\r
 353   gFile->Close();\r
 354   return kTRUE;\r
 355 }\r
 356 \r
 357 //________________________________________________________\r
 358 Bool_t AliTRDrecoTask::LoadDetectorMap(const Char_t *file, const Char_t *dir)\r
 359 {\r
 360 // Load detector map.\r
 361 \r
 362   if(!TFile::Open(file)){\r
 363     AliWarning(Form("Couldn't open file %s.", file));\r
 364     return kFALSE;\r
 365   }\r
 366   if(!gFile->cd(dir)){\r
 367     AliWarning(Form("Couldn't cd to %s in %s.", dir, file));\r
 368     gFile->Close();\r
 369     return kFALSE;\r
 370   }\r
 371   TObjArray *info = NULL;\r
 372   if(!(info = (TObjArray*)gDirectory->Get("TRDinfoGen"))){\r
 373     AliWarning("Missing TRDinfoGen container.");\r
 374     gFile->Close();\r
 375     return kFALSE;\r
 376   }\r
 377   fDets = (TObjArray*)((TObjArray*)info->FindObject("Chambers Status"))->Clone();\r
 378   if(!fDets){\r
 379     if(!info->At(4) || strcmp("TObjArray", info->At(4)->IsA()->GetName())) AliError("Looking for old style chamber status map. Failed.");\r
 380     else {\r
 381       AliWarning("Looking for old style chamber status map.");\r
 382       fDetsV = (TObjArray*)((TObjArray*)info->At(4))->Clone();\r
 383     }\r
 384   }\r
 385   gFile->Close();\r
 386   info->Delete(); delete info;\r
 387   return kTRUE;\r
 388 }\r
 389 \r
 390 \r
 391 //________________________________________________________\r
 392 Bool_t AliTRDrecoTask::Save(TObjArray * const results){\r
 393   //\r
 394   // Store the output graphs in a ROOT file\r
 395   // Input TObject array will not be written as Key to the file,\r
 396   // only content itself\r
 397   //\r
 398 \r
 399   TDirectory *cwd = gDirectory;\r
 400   if(!TFile::Open(Form("TRD.Result%s.root", GetName()), "RECREATE")) return kFALSE;\r
 401 \r
 402   TIterator *iter = results->MakeIterator();\r
 403   TObject *inObject = NULL, *outObject = NULL;\r
 404   while((inObject = iter->Next())){\r
 405     outObject = inObject->Clone();\r
 406     outObject->Write(NULL, TObject::kSingleKey);\r
 407   }\r
 408   delete iter;\r
 409   gFile->Close(); delete gFile;\r
 410   cwd->cd(); \r
 411   return kTRUE;\r
 412 }\r
 413 \r
 414 //_______________________________________________________\r
 415 Bool_t AliTRDrecoTask::PostProcess()\r
 416 {\r
 417 // To be implemented by particular tasks\r
 418 \r
 419   AliWarning("Post processing of reference histograms not implemented.");\r
 420   return kTRUE;\r
 421 }\r
 422 \r
 423 //_______________________________________________________\r
 424 void AliTRDrecoTask::MakeDetectorPlot(Int_t ly, const Option_t *opt)\r
 425 {\r
 426 // Draw chamber boundaries in eta/phi plots with misalignments\r
 427 // based on info collected by AliTRDinfoGen\r
 428 \r
 429   if(!fDets){\r
 430     AliWarning("NEW Detector map and status not available. Try OLD");\r
 431     MakeDetectorPlotOLD(ly, opt);\r
 432     return;\r
 433   }\r
 434   AliTRDchmbInfo *ci(NULL);\r
 435   for(Int_t idet(0); idet<fDets->GetEntriesFast(); idet++){\r
 436     if(!(ci = (AliTRDchmbInfo*)fDets->At(idet))) continue;\r
 437     if(AliTRDgeometry::GetLayer(ci->GetDetector()) != ly) continue;\r
 438     ci->Draw(opt);\r
 439   }\r
 440 }\r
 441 \r
 442 //_______________________________________________________\r
 443 void AliTRDrecoTask::MakeDetectorPlotOLD(Int_t ly, const Option_t *opt)\r
 444 {\r
 445 // Draw chamber boundaries in eta/phi plots with misalignments\r
 446 // based on info collected by AliTRDinfoGen OLD data storage\r
 447 \r
 448   if(!fDetsV){\r
 449     AliError("OLD Detector map and status not available.");\r
 450     return;\r
 451   }\r
 452   if(!fDetsV->GetEntries()){\r
 453     AliError("OLD Detector map and status not filled.");\r
 454     return;\r
 455   }\r
 456 \r
 457   Float_t xmin(0.), xmax(0.);\r
 458   TBox *gdet = new TBox();\r
 459   gdet->SetLineColor(kBlack);gdet->SetFillColor(kBlack);\r
 460   Int_t style[] = {0, 3003};\r
 461   for(Int_t idet(0); idet<540; idet++){\r
 462     if(idet%6 != ly) continue;\r
 463     TVectorF *det((TVectorF*)fDetsV->At(idet));\r
 464     if(!det) continue;\r
 465     Int_t iopt = Int_t((*det)[4]);\r
 466     if(strcmp(opt, "eta")==0){\r
 467       xmin=(*det)[0]; xmax=(*det)[2];\r
 468     } else if(strcmp(opt, "pad")==0){\r
 469       Int_t stk(AliTRDgeometry::GetStack(idet));\r
 470       xmin=-0.6+16*(4-stk)-(stk<2?4:0); xmax=xmin+(stk==2?12:16)-0.2;\r
 471     } else continue;\r
 472     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
 473     if(iopt==1){\r
 474       gdet->SetFillStyle(style[1]);gdet->SetFillColor(kBlack);\r
 475       gdet->DrawBox(xmin, (*det)[1], xmax, (*det)[3]);\r
 476     } else {\r
 477       gdet->SetFillStyle(style[0]);\r
 478       gdet->DrawBox(xmin, (*det)[1], xmax, (*det)[3]);\r
 479       if(iopt==2){\r
 480         gdet->SetFillStyle(style[1]);gdet->SetFillColor(kGreen);\r
 481         gdet->DrawBox(xmin, (*det)[1], xmax, 0.5*((*det)[3]+(*det)[1]));\r
 482       } else if(iopt==3){\r
 483         gdet->SetFillStyle(style[1]);gdet->SetFillColor(kRed);\r
 484         gdet->DrawBox(xmin, 0.5*((*det)[3]+(*det)[1]), xmax, (*det)[3]);\r
 485       } else if(iopt!=0) AliError(Form("Wrong chmb. status[%d] for det[%03d]", iopt, idet));\r
 486     }\r
 487   }\r
 488   Float_t dsm = TMath::TwoPi()/AliTRDgeometry::kNsector;\r
 489   xmin=0.;\r
 490   if(strcmp(opt, "pad")==0) xmin=38.;\r
 491   TLatex *sm = new TLatex(); sm->SetTextAlign(22);sm->SetTextColor(kBlack); sm->SetTextFont(32);sm->SetTextSize(0.03);\r
 492   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
 493 }\r
 494 \r
 495 \r
 496 //_______________________________________________________\r
 497 void AliTRDrecoTask::MakeSummary()\r
 498 {\r
 499 // To be implemented by particular tasks\r
 500   AliWarning("Summary not available");\r
 501 }\r
 502 \r
 503 //_______________________________________________________\r
 504 void AliTRDrecoTask::SetDebugLevel(Int_t level)\r
 505 {\r
 506 // Generic debug handler\r
 507 \r
 508   AliAnalysisTaskSE::SetDebugLevel(level);\r
 509   if(DebugLevel()>=1){\r
 510     AliInfo(Form("Debug Level for Task %s set to %d", GetName(), level));\r
 511     TDirectory *savedir = gDirectory;\r
 512     fgDebugStream = new TTreeSRedirector("TRD.DebugPerformance.root");\r
 513     savedir->cd();\r
 514   }\r
 515 }\r
 516 \r
 517 //____________________________________________________________________\r
 518 void AliTRDrecoTask::Terminate(Option_t *)\r
 519 {\r
 520   //\r
 521   // Terminate\r
 522   //\r
 523 \r
 524   if(fgDebugStream){ \r
 525     delete fgDebugStream;\r
 526     fgDebugStream = NULL;\r
 527   }\r
 528   fContainer = dynamic_cast<TObjArray *>(GetOutputData(1));\r
 529   if(fContainer && fRunTerminate){\r
 530     PostProcess();\r
 531     MakeSummary();\r
 532   }\r
 533 }\r
 534 \r
 535 //________________________________________________________\r
 536 Float_t AliTRDrecoTask::SetNormZ(TH2 *h2, Int_t bxmin, Int_t bxmax, Int_t bymin, Int_t bymax, Float_t thr)\r
 537 {\r
 538 // Normalize histo content to the mean value in the range specified by bin ranges\r
 539 // [bxmin, bxmax] on the x axis and [bymin, bymax] on the y axis.\r
 540 // Optionally a threshold "thr" can be specified to disregard entries with no meaning\r
 541 \r
 542   Float_t s = 0., c=0.; Int_t is(0);\r
 543   for(Int_t ix(bxmin); ix<=(bxmax>0?bxmax:(h2->GetXaxis()->GetNbins())); ix++){\r
 544     for(Int_t iy(bymin); iy<=(bymax>0?bymax:(h2->GetYaxis()->GetNbins())); iy++){\r
 545       if((c = h2->GetBinContent(ix, iy))<thr) continue;\r
 546       s += c; is++;\r
 547     }\r
 548   }\r
 549   s/= (is?is:1);\r
 550   for(Int_t ix(1); ix<=h2->GetXaxis()->GetNbins(); ix++){\r
 551     for(Int_t iy(1); iy<=h2->GetYaxis()->GetNbins(); iy++){\r
 552       if((c = h2->GetBinContent(ix, iy))<thr) h2->SetBinContent(ix, iy, thr-1000);\r
 553       else h2->SetBinContent(ix, iy, 100.*(c/s-1.));\r
 554     }\r
 555   }\r
 556   return s;\r
 557 }\r
 558 \r
 559 //________________________________________________________\r
 560 void AliTRDrecoTask::SetRangeZ(TH2 *h2, Float_t min, Float_t max, Float_t thr)\r
 561 {\r
 562 // Set range on Z axis such to avoid outliers\r
 563 \r
 564   Float_t c(0.), dz(1.e-3*(max-min));\r
 565   for(Int_t ix(1); ix<=h2->GetXaxis()->GetNbins(); ix++){\r
 566     for(Int_t iy(1); iy<=h2->GetYaxis()->GetNbins(); iy++){\r
 567       if((c = h2->GetBinContent(ix, iy))<thr) continue;\r
 568       if(c<=min) h2->SetBinContent(ix, iy, min+dz);\r
 569     }\r
 570   }\r
 571   h2->GetZaxis()->SetRangeUser(min, max);\r
 572 }\r
 573 \r
 574 //________________________________________________________\r
 575 Float_t AliTRDrecoTask::GetMeanStat(TH1 *h, Float_t cut, Int_t opt, Float_t *sigma)\r
 576 {\r
 577 // Return mean number of entries/bin of histogram "h".\r
 578 // If optionally sigma is allocated than it is also filled with sigma paramter of the gauss fit \r
 579 //\r
 580 // Option "opt" is given the following values are accepted:\r
 581 //   -1 : consider only entries less than "cut"\r
 582 //   1  : consider only entries greater than "cut"\r
 583 //   0  : no "cut" [dafault]\r
 584 // Error codes\r
 585 //   -999. : statistics too low [20]\r
 586 //   -998. : fit failed\r
 587 \r
 588   const Int_t kvd(100000);\r
 589   Float_t v[kvd];\r
 590   Int_t nbx(h->GetNbinsX()), nby(h->GetNbinsY()), nbz(h->GetNbinsZ());\r
 591   Int_t nv(0); Float_t xmin(1.e5), xmax(-xmin);\r
 592   for(Int_t ix(1); ix<=nbx; ix++){\r
 593     for(Int_t iy(1); iy<=nby; iy++){\r
 594       for(Int_t iz(1); iz<=nbz; iz++){\r
 595         Float_t c = h->GetBinContent(ix, iy, iz);\r
 596         if(opt*(c-cut)<0.) continue;\r
 597         v[nv++] = c;\r
 598         if(c<xmin) xmin = c;\r
 599         if(c>xmax) xmax = c;\r
 600         if(nv==kvd){\r
 601           printf("W - AliTRDrecoTask::GetMeanStat() :: Unreliable results for %s[%s]. Statical allocation exceeded.\n", h->GetName(), h->GetTitle());\r
 602           break;\r
 603         }\r
 604       }\r
 605       if(nv==kvd) break;\r
 606     }\r
 607     if(nv==kvd) break;\r
 608   }\r
 609   if(nv<10){\r
 610     //printf("W - AliTRDrecoTask::GetMeanStat() :: Failed for %s[%s]. Statical undefined [%d].\n", h->GetName(), h->GetTitle(), nv);\r
 611     return -999.;\r
 612   }\r
 613   if(fgProjector) delete fgProjector;\r
 614   fgProjector = new TH1F("hProjector", "", 20, 0.5*(3*xmin-xmax), 0.5*(3*xmax - xmin));\r
 615   for(Int_t iv(0); iv<nv; iv++) fgProjector->Fill(v[iv]);\r
 616   TF1 f("f", "gaus", xmin, xmax);\r
 617   f.SetParameter(0, fgProjector->Integral());\r
 618   f.SetParameter(1, fgProjector->GetMean()); f.SetParLimits(1, xmin, xmax);\r
 619   f.SetParameter(2, fgProjector->GetRMS());\r
 620   if(fgProjector->Fit(&f, "WQ0", "goff")) return -998.;\r
 621   if(sigma) *sigma = f.GetParameter(2);\r
 622   return f.GetParameter(1);\r
 623 }\r
 624 \r
 625 //________________________________________________________\r
 626 Int_t AliTRDrecoTask::Rebin(TH2 *h, Int_t n, Int_t rebinX[], Int_t rebinY[], Int_t nstat)\r
 627 {\r
 628 // Rebin histo "h" according to "rebinning" strategy "n" steps such to obtain mean statistics per bin over "nstat"\r
 629 \r
 630   Int_t irebin(0);\r
 631   while(irebin<n && GetMeanStat(h, .5, 1)<nstat){\r
 632     h->Rebin2D(rebinX[irebin], rebinY[irebin]);\r
 633     irebin++;\r
 634   }\r
 635   return irebin;\r
 636 }\r
 637 \r
 638 //________________________________________________________\r
 639 AliTRDrecoTask::AliTRDrecoProjection::AliTRDrecoProjection()\r
 640   :TNamed()\r
 641   ,fH(NULL)\r
 642   ,fNrebin(0)\r
 643 {\r
 644   // constructor\r
 645   fRebin[0] = NULL;fRebin[1] = NULL;\r
 646   memset(fAx, 0, 3*sizeof(Int_t));\r
 647   memset(fRange, 0, 4*sizeof(Float_t));\r
 648 }\r
 649 \r
 650 //________________________________________________________\r
 651 AliTRDrecoTask::AliTRDrecoProjection::~AliTRDrecoProjection()\r
 652 {\r
 653   // destructor\r
 654   if(fH) delete fH;\r
 655   if(fRebin[0]) delete [] fRebin[0];\r
 656   if(fRebin[1]) delete [] fRebin[1];\r
 657 }\r
 658 \r
 659 //________________________________________________________\r
 660 void AliTRDrecoTask::AliTRDrecoProjection::Build(const Char_t *n, const Char_t *t, Int_t ix, Int_t iy, Int_t iz, TAxis *aa[])\r
 661 {\r
 662 // check and build (if neccessary) projection determined by axis "ix", "iy" and "iz"\r
 663   if(!aa[ix] || !aa[iy] || !aa[iz]) return;\r
 664   TAxis *ax(aa[ix]), *ay(aa[iy]), *az(aa[iz]);\r
 665   // check ax definiton to protect against older versions of the data\r
 666   if(ax->GetNbins()<=0 || (ax->GetXmax()-ax->GetXmin())<=0.){\r
 667     AliWarning(Form("Wrong definition of axis[%d] \"%s\"[%d](%f %f).", ix, ax->GetTitle(), ax->GetNbins(), ax->GetXmin(), ax->GetXmax()));\r
 668     return;\r
 669   }\r
 670   if(ay->GetNbins()<=0 || (ay->GetXmax()-ay->GetXmin())<=0.){\r
 671     AliWarning(Form("Wrong definition of axis[%d] \"%s\"[%d](%f %f).", ix, ay->GetTitle(), ay->GetNbins(), ay->GetXmin(), ay->GetXmax()));\r
 672     return;\r
 673   }\r
 674   if(az->GetNbins()<=0 || (az->GetXmax()-az->GetXmin())<=0.){\r
 675     AliWarning(Form("Wrong definition of axis[%d] \"%s\"[%d](%f %f).", ix, az->GetTitle(), az->GetNbins(), az->GetXmin(), az->GetXmax()));\r
 676     return;\r
 677   }\r
 678   SetNameTitle(n,t);\r
 679   fH = new TH3I(n, Form("%s;%s;%s;%s", t, ax->GetTitle(), ay->GetTitle(), az->GetTitle()),\r
 680     ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),\r
 681     ay->GetNbins(), ay->GetXmin(), ay->GetXmax(),\r
 682     az->GetNbins(), az->GetXmin(), az->GetXmax());\r
 683   fAx[0] = ix; fAx[1] = iy; fAx[2] = iz;\r
 684   fRange[0] = az->GetXmin()/3.; fRange[1] = az->GetXmax()/3.;\r
 685   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
 686     ax->GetTitle(), ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),\r
 687     ay->GetTitle(), ay->GetNbins(), ay->GetXmin(), ay->GetXmax(),\r
 688     az->GetTitle(), az->GetNbins(), az->GetXmin(), az->GetXmax()));\r
 689 }\r
 690 \r
 691 //________________________________________________________\r
 692 AliTRDrecoTask::AliTRDrecoProjection& AliTRDrecoTask::AliTRDrecoProjection::operator=(const AliTRDrecoProjection& rhs)\r
 693 {\r
 694 // copy projections\r
 695   if(this == &rhs) return *this;\r
 696 \r
 697   TNamed::operator=(rhs);\r
 698   if(fNrebin){fNrebin=0; delete [] fRebin[0]; delete [] fRebin[1];}\r
 699   if(rhs.fNrebin) SetRebinStrategy(rhs.fNrebin, rhs.fRebin[0], rhs.fRebin[1]);\r
 700   memcpy(fAx, rhs.fAx, 3*sizeof(Int_t));\r
 701   memcpy(fRange, rhs.fRange, 4*sizeof(Float_t));\r
 702   if(fH) delete fH;\r
 703   if(rhs.fH) fH=(TH3I*)rhs.fH->Clone(Form("%s_CLONE", rhs.fH->GetName()));\r
 704   return *this;\r
 705 }\r
 706 \r
 707 //________________________________________________________\r
 708 AliTRDrecoTask::AliTRDrecoProjection& AliTRDrecoTask::AliTRDrecoProjection::operator+=(const AliTRDrecoProjection& other)\r
 709 {\r
 710 // increment projections\r
 711   if(!fH || !other.fH) return *this;\r
 712   AliDebug(2, Form("%s+=%s [%s+=%s]", GetName(), other.GetName(), fH->GetName(), (other.fH)->GetName()));\r
 713   fH->Add(other.fH);\r
 714   return *this;\r
 715 }\r
 716 \r
 717 //________________________________________________________\r
 718 void AliTRDrecoTask::AliTRDrecoProjection::Increment(Int_t bin[], Double_t v)\r
 719 {\r
 720 // increment bin with value "v" pointed by general coord in "bin"\r
 721   if(!fH) return;\r
 722   AliDebug(4, Form("  %s[%2d]", fH->GetName(), Int_t(v)));\r
 723   fH->AddBinContent(fH->GetBin(bin[fAx[0]],bin[fAx[1]],bin[fAx[2]]), Int_t(v));\r
 724 }\r
 725 \r
 726 //________________________________________________________\r
 727 Double_t AliTRDrecoTask::AliTRDrecoProjection::GetTrendValue(const Int_t mid, Double_t *e, Double_t *s, Double_t *se) const\r
 728 {\r
 729 //   Return result of fitting the main distribution (represented on the z axis) with the function selected\r
 730 // "mid". Optionally return the Mean and RMS of the distribution pointing to "m" and "s"\r
 731 \r
 732   if(!fH){\r
 733     AliDebug(1, Form("Missing 3D in %s", GetName()));\r
 734     return -999.;\r
 735   }\r
 736   TH1 *h1s(NULL);\r
 737   if(!(h1s = (TH1D*)fH->Project3D("z"))){\r
 738     AliDebug(1, Form("Failed Project3D(\"z\") in %s", GetName()));\r
 739     return -999.;\r
 740   }\r
 741   Int_t ne((Int_t)h1s->Integral());\r
 742   if(ne<30){\r
 743     AliDebug(1, Form("Statistics too low[%2d] in %s", ne, GetName()));\r
 744     delete h1s;\r
 745     return -999.;\r
 746   }\r
 747   TAxis *az(h1s->GetXaxis());\r
 748   Float_t mn(h1s->GetMean()), rms(h1s->GetRMS()),\r
 749           v(mn),  // main trending value (mean, mu, MPV)\r
 750           ve(rms),// dispersion (RMS, sigma, landau 2nd param)\r
 751           ev(0.), // error on v\r
 752           eve(0.);// error on ve\r
 753   if(mid==1){\r
 754     TF1 fg("fg", "gaus", az->GetXmin(), az->GetXmax());\r
 755     fg.SetParameter(0, Float_t(ne)); fg.SetParameter(1, mn); fg.SetParameter(2, rms);\r
 756     h1s->Fit(&fg, "WQ0");\r
 757     v = fg.GetParameter(1); ev = fg.GetParError(1);\r
 758     ve= fg.GetParameter(2);\r
 759   } else if (mid==2) {\r
 760     TF1 fl("fl", "landau", az->GetXmin(), az->GetXmax());\r
 761     fl.SetParameter(0, Float_t(ne)); fl.SetParameter(1, mn); fl.SetParameter(2, rms);\r
 762     h1s->Fit(&fl, "WQ0");\r
 763     v = fl.GetMaximumX(); ev = fl.GetParError(1);\r
 764     ve= fl.GetParameter(2);\r
 765   }\r
 766   if(e)  *e  = ev;\r
 767   if(s)  *s  = ve;\r
 768   if(se) *se = eve;\r
 769   AliDebug(2, Form("[%d] %s(%4d) = M{%f+-%f} S{%f+-%f}", mid, fH->GetName(), (Int_t)h1s->Integral(), v, ev, ve, eve));\r
 770 \r
 771   delete h1s;\r
 772   return v;\r
 773 }\r
 774 \r
 775 //________________________________________________________\r
 776 TH2* AliTRDrecoTask::AliTRDrecoProjection::Projection2Dbin(Int_t bin, Bool_t mc)\r
 777 {\r
 778 // dumb 2D projection for bin including under/over flow. Default all [bin==-1]\r
 779 \r
 780   TAxis *ax(fH->GetXaxis()), *ay(fH->GetYaxis()), *az(fH->GetZaxis());\r
 781   Int_t nbins(az->GetNbins());\r
 782   TH2F *h2(NULL);\r
 783   if(bin<0) h2 = new TH2F(Form("%s_2D", fH->GetName()),\r
 784                 Form("%s;%s;%s;Entries", fH->GetTitle(), ax->GetTitle(), ay->GetTitle()),\r
 785                 ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),\r
 786                 ay->GetNbins(), ay->GetXmin(), ay->GetXmax());\r
 787   else h2 = new TH2F(Form("%s%d_2D", fH->GetName(), bin),\r
 788                 Form("%s | #it{%4.2f<=p_{t}^{%s}[GeV/c]<%4.2f};%s;%s;Entries", fH->GetTitle(),\r
 789                 bin?fgPt[bin-1]:0., mc?"MC":"", bin>nbins?99.99:fgPt[bin], ax->GetTitle(), ay->GetTitle()),\r
 790                 ax->GetNbins(), ax->GetXmin(), ax->GetXmax(),\r
 791                 ay->GetNbins(), ay->GetXmin(), ay->GetXmax());\r
 792   for(Int_t ix(1); ix<=ax->GetNbins(); ix++){\r
 793     for(Int_t iy(1); iy<=ay->GetNbins(); iy++){\r
 794       Int_t ibin = h2->GetBin(ix, iy);\r
 795       for(Int_t iz(0); iz<=az->GetNbins()+1; iz++){\r
 796         if(bin<0) h2->AddBinContent(ibin, fH->GetBinContent(ix, iy, iz));\r
 797         else if(bin==iz){\r
 798           h2->AddBinContent(ibin, fH->GetBinContent(ix, iy, iz));\r
 799           break;\r
 800         }\r
 801       }\r
 802     }\r
 803   }\r
 804   return h2;\r
 805 }\r
 806 \r
 807 //________________________________________________________\r
 808 TH2* AliTRDrecoTask::AliTRDrecoProjection::Projection2D(const Int_t nstat, const Int_t ncol, const Int_t mid, Bool_t del)\r
 809 {\r
 810 // build the 2D projection and adjust binning\r
 811 \r
 812   const Char_t *title[] = {"Mean", "#mu", "MPV"};\r
 813   if(!fH){\r
 814     AliDebug(1, Form("Missing 3D in %s", GetName()));\r
 815     return NULL;\r
 816   }\r
 817   TAxis *ax(fH->GetXaxis()), *ay(fH->GetYaxis()), *az(fH->GetZaxis());\r
 818   TH2D *h2s(NULL), *hyx(NULL);\r
 819   if(!(h2s = (TH2D*)fH->Project3D("yx"))){\r
 820     AliDebug(1, Form("Failed Project3D(\"yx\") in %s", GetName()));\r
 821     return NULL;\r
 822   }\r
 823   // save a copy of the original distribution\r
 824   if(!del){\r
 825     hyx = (TH2D*)h2s->Clone();\r
 826     hyx->SetName(Form("%sEn", fH->GetName()));\r
 827   }\r
 828   Int_t irebin(Rebin(h2s, fNrebin, fRebin[0], fRebin[1], nstat)), dxBin(1), dyBin(1);\r
 829   for(Int_t ir(0); ir<irebin; ir++){dxBin*=fRebin[0][ir]; dyBin*=fRebin[1][ir];}\r
 830   Int_t nx(h2s->GetNbinsX()), ny(h2s->GetNbinsY());\r
 831   delete h2s;\r
 832   if(mid<0) return NULL;\r
 833 \r
 834   // start projection\r
 835   TH1 *h(NULL); Int_t n(0);\r
 836   Float_t dz=(fRange[1]-fRange[1])/ncol;\r
 837   TString titlez(az->GetTitle()); TObjArray *tokenTitle(titlez.Tokenize(" "));\r
 838   Int_t nt(tokenTitle->GetEntriesFast());\r
 839   TH2 *h2(NULL);\r
 840   if((h2 = (TH2*)gDirectory->Get(Form("%s_2D", fH->GetName())))) delete h2; // avoid ROOT warning messages\r
 841   h2 = new TH2F(Form("%s_2D", fH->GetName()),\r
 842             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
 843             nx, ax->GetXmin(), ax->GetXmax(), ny, ay->GetXmin(), ay->GetXmax());\r
 844   tokenTitle->Delete(); delete tokenTitle;\r
 845   h2->SetContour(ncol);\r
 846   h2->GetZaxis()->CenterTitle();\r
 847   h2->GetZaxis()->SetTitleOffset(1.4);\r
 848   h2->GetZaxis()->SetRangeUser(fRange[0], fRange[1]);\r
 849   AliDebug(2, Form("%s[%s] nx[%d] ny[%d]", h2->GetName(), h2->GetTitle(), nx, ny));\r
 850   for(Int_t iy(0); iy<ny; iy++){\r
 851     for(Int_t ix(0); ix<nx; ix++){\r
 852       h = fH->ProjectionZ(Form("%s_z", h2->GetName()), ix*dxBin+1, (ix+1)*dxBin, iy*dyBin+1, (iy+1)*dyBin);\r
 853       Int_t ne((Int_t)h->Integral());\r
 854       //printf("  x[%2d %2d] y[%2d %2d] ne[%4d]\n", ix*dxBin+1, (ix+1)*dxBin, iy*dyBin+1, (iy+1)*dyBin, ne);\r
 855       if(ne<nstat/2){\r
 856         h2->SetBinContent(ix+1, iy+1, -999);\r
 857         h2->SetBinError(ix+1, iy+1, 1.);\r
 858         n++;\r
 859       }else{\r
 860         // redo the projection by adding 1 bin @ left and 1 bin @ right for smoothing\r
 861         h = fH->ProjectionZ(Form("%s_z", h2->GetName()), ix*dxBin, (ix+1)*dxBin+1, iy*dyBin, (iy+1)*dyBin+1);\r
 862         Float_t v(h->GetMean()), ve(h->GetRMS());\r
 863         if(mid==1){\r
 864           TF1 fg("fg", "gaus", az->GetXmin(), az->GetXmax());\r
 865           fg.SetParameter(0, Float_t(ne)); fg.SetParameter(1, v); fg.SetParameter(2, ve);\r
 866           h->Fit(&fg, "WQ0");\r
 867           v = fg.GetParameter(1); ve = fg.GetParameter(2);\r
 868         } else if (mid==2) {\r
 869           TF1 fl("fl", "landau", az->GetXmin(), az->GetXmax());\r
 870           fl.SetParameter(0, Float_t(ne)); fl.SetParameter(1, v); fl.SetParameter(2, ve);\r
 871           h->Fit(&fl, "WQ0");\r
 872           v = fl.GetMaximumX(); ve = fl.GetParameter(2);\r
 873 /*          TF1 fgle("gle", "[0]*TMath::Landau(x, [1], [2], 1)*TMath::Exp(-[3]*x/[1])", az->GetXmin(), az->GetXmax());\r
 874           fgle.SetParameter(0, fl.GetParameter(0));\r
 875           fgle.SetParameter(1, fl.GetParameter(1));\r
 876           fgle.SetParameter(2, fl.GetParameter(2));\r
 877           fgle.SetParameter(3, 1.);fgle.SetParLimits(3, 0., 5.);\r
 878           h->Fit(&fgle, "WQ");\r
 879           v = fgle.GetMaximumX(); ve = fgle.GetParameter(2);*/\r
 880         }\r
 881         if(v<fRange[0]) h2->SetBinContent(ix+1, iy+1, fRange[0]+0.1*dz);\r
 882         else h2->SetBinContent(ix+1, iy+1, v);\r
 883         h2->SetBinError(ix+1, iy+1, ve);\r
 884       }\r
 885     }\r
 886   }\r
 887   if(h) delete h;\r
 888   if(n==nx*ny){  // clean empty projections\r
 889     AliDebug(1, Form("Empty projection in %s", GetName()));\r
 890     delete h2; h2=NULL;\r
 891   }\r
 892   return h2;\r
 893 }\r
 894 \r
 895 //________________________________________________________\r
 896 void AliTRDrecoTask::AliTRDrecoProjection::SetRebinStrategy(Int_t n, Int_t rebx[], Int_t reby[])\r
 897 {\r
 898 // define rebinning strategy for this projection\r
 899   fNrebin = n;\r
 900   fRebin[0] = new Int_t[n]; memcpy(fRebin[0], rebx, n*sizeof(Int_t));\r
 901   fRebin[1] = new Int_t[n]; memcpy(fRebin[1], reby, n*sizeof(Int_t));\r
 902 }\r
 903 \r
 904 //________________________________________________________\r
 905 void AliTRDrecoTask::SetTriggerList(const Char_t *tl)\r
 906 {\r
 907 // Store list of triggers to be monitored\r
 908   TString stl(tl);\r
 909   if(fTriggerList){ fTriggerList->Delete(); delete fTriggerList;}\r
 910   TObjArray *atl = stl.Tokenize(" ");\r
 911   fTriggerList = (TObjArray*)atl->Clone("");\r
 912   atl->Delete(); delete atl;\r
 913   AliInfo("Running only for triggers::");\r
 914   fTriggerList->Print();\r
 915 }\r
 916 \r
 917 \r