PWGPP/TRD/AliTRDrecoTask.cxx

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