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