PWGUD/multVScentPbPb/AliTrackletTaskMulti.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, ALICE Experiment at CERN, All rights reserved. *
   3 *                                                                        *
   4 * Author: The ALICE Off-line Project.                                    *
   5 * Contributors are mentioned in the code where appropriate.              *
   6 *                                                                        *
   7 * Permission to use, copy, modify and distribute this software and its   *
   8 * documentation strictly for non-commercial purposes is hereby granted   *
   9 * without fee, provided that the above copyright notice appears in all   *
  10 * copies and that both the copyright notice and this permission notice   *
  11 * appear in the supporting documentation. The authors make no claims     *
  12 * about the suitability of this software for any purpose. It is          *
  13 * provided "as is" without express or implied warranty.                  *
  14 **************************************************************************/
  15
  16 ///////////////////////////////////////////////////////////////////////////
  17 // Class AliTrackletTaskMulti                                            //
  18 // Analysis task to produce data and MC histos needed for tracklets      //
  19 // dNdEta extraction in multiple bins in one go                          //
  20 // Author:  ruben.shahoyan@cern.ch                                       //
  21 ///////////////////////////////////////////////////////////////////////////
  22 /*
  23   Important parameters to set:
  24   1) make sure to initialize correct geometry in UserCreateOutputObjects
  25   2) The cut on signal selection variable (delta, dphi ...) should be decided beforehand
  26 ...
  27 */
  28
  29 #include "TFile.h"
  30 #include "TChain.h"
  31 #include "TTree.h"
  32 #include "TRandom.h"
  33 #include "TH1F.h"
  34 #include "TH2F.h"
  35 #include "TH3F.h"
  36 #include "THnSparse.h"
  37 #include "TList.h"
  38 #include "TNtuple.h"
  39 #include "TObjArray.h"
  40 #include "TGeoGlobalMagField.h"
  41
  42 #include "AliAnalysisManager.h"
  43
  44 #include "AliMultiplicity.h"
  45 #include "AliESDEvent.h"
  46 #include "AliESDInputHandler.h"
  47 #include "AliESDInputHandlerRP.h"
  48 #include "../ANALYSIS/EventMixing/AliMixEventInputHandler.h"
  49 #include "AliCDBPath.h"
  50 #include "AliCDBManager.h"
  51 #include "AliCDBEntry.h"
  52 #include "AliCDBStorage.h"
  53 #include "AliGeomManager.h"
  54 #include "AliMagF.h"
  55 #include "AliESDVZERO.h"
  56 #include "AliESDZDC.h"
  57 #include "AliRunLoader.h"
  58 #include "AliMCEventHandler.h"
  59 #include "AliMCEvent.h"
  60 #include "AliMCParticle.h"
  61 #include "AliStack.h"
  62 #include "AliGenEventHeader.h"
  63 #include "AliCentrality.h"
  64 #include "../ITS/AliITSRecPoint.h"
  65 #include "../ITS/AliITSgeomTGeo.h"
  66 #include "../ITS/AliITSMultReconstructor.h"
  67
  68 #include "AliLog.h"
  69
  70 #include "AliPhysicsSelection.h"
  71 #include "AliTrackletTaskMulti.h"
  72 #include "AliITSMultRecBg.h"
  73 #include "AliGenEventHeader.h"
  74 #include "AliGenHijingEventHeader.h"
  75 #include "AliGenDPMjetEventHeader.h"
  76 #include "AliESDtrackCuts.h"
  77
  78 ClassImp(AliTrackletTaskMulti)
  79
  80 // centrality percentile (inverted: 100% - most central)
  81 //const Float_t  AliTrackletTaskMulti::fgkCentPerc[] = {0,100};//{0,5,10,20,30};
  82 const Float_t  AliTrackletTaskMulti::fgkCentPerc[] = {0,5,10,20,30,40,50,60,70,80,90};
  83   //{0,10,20,30,40,50,60,70,80,90,95,101};
  84
  85 const char* AliTrackletTaskMulti::fgCentSelName[] = {"V0M","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","TKLvsV0M","ZEMvsZDC"};
  86
  87 const char*  AliTrackletTaskMulti::fgkPDGNames[] = {
  88 "#pi^{+}",
  89 "p",
  90 "K^{+}",
  91 "K^{*+}",
  92 "e^{-}",
  93 "#mu^{-}",
  94 "#rho^{+}",
  95 "D^{+}",
  96 "D^{*+}",
  97 "D_{s}^{+}",
  98 "D_{s}^{*+}",
  99 "#Delta^{-}",
 100 "#Delta^{+}",
 101 "#Delta^{++}",
 102 "#Sigma^{-}",
 103 "#Sigma^{+}",
 104 "#Sigma^{*-}",
 105 "#Sigma^{*+}",
 106 "#Sigma^{*+}_{c}",
 107 "#Sigma^{*++}_{c}",
 108 "#Xi^{-}",
 109 "#Xi^{*-}",
 110 "#Lambda^{+}_{c}",
 111 "n",
 112 "#Delta^{0}",
 113 "#gamma",
 114 "K^{0}_{S}",
 115 "K^{0}_{L}",
 116 "K^{0}",
 117 "K^{*}",
 118 "#eta",
 119 "#pi^{0}",
 120 "#rho^{0}",
 121 "#varphi",
 122 "#eta'",
 123 "#omega",
 124 "#Lambda",
 125 "#Sigma^{0}",
 126 "#Sigma^{*0}_{c}",
 127 "#Sigma^{*0}",
 128 "D^{0}",
 129 "D^{*0}",
 130 "#Xi_{0}",
 131 "#Xi^{*0}",
 132 "#Xi^{0}_{c}",
 133 "#Xi^{*0}_{c}",
 134 "Nuclei",
 135 "Others"
 136 };
 137
 138 const int AliTrackletTaskMulti::fgkPDGCodes[] = {
 139   211,
 140  2212,
 141   321,
 142   323,
 143    11,
 144    13,
 145   213,
 146   411,
 147   413,
 148   431,
 149   433,
 150  1114,
 151  2214,
 152  2224,
 153  3112,
 154  3222,
 155  3114,
 156  3224,
 157  4214,
 158  4224,
 159  3312,
 160  3314,
 161  4122,
 162  2112,
 163  2114,
 164    22,
 165   310,
 166   130,
 167   311,
 168   313,
 169   221,
 170   111,
 171   113,
 172   333,
 173   331,
 174   223,
 175  3122,
 176  3212,
 177  4114,
 178  3214,
 179   421,
 180   423,
 181  3322,
 182  3324,
 183  4132,
 184  4314
 185 // nuclei
 186 // unknown
 187 };
 188
 189 const double kEtaBinWidth=0.25;
 190
 191 //________________________________________________________________________
 192 /*//Default constructor
 193 AliTrackletTaskMulti::AliTrackletTaskMulti(const char *name)
 194   : AliAnalysisTaskSE(name),
 195 */
 196 //________________________________________________________________________
 197 AliTrackletTaskMulti::AliTrackletTaskMulti(const char *name)
 198   : AliAnalysisTaskSE(name),
 199 //
 200   fOutput(0),
 201 //
 202   fDoNormalReco(kFALSE),
 203   fDoInjection(kFALSE),
 204   fDoRotation(kFALSE),
 205   fDoMixing(kFALSE),
 206   //
 207   fUseMC(kFALSE),
 208   fCheckReconstructables(kFALSE),
 209 //
 210   fHistosTrData(0),
 211   fHistosTrInj(0),
 212   fHistosTrRot(0),
 213   fHistosTrMix(0),
 214 //
 215   fHistosTrPrim(0),
 216   fHistosTrSec(0),
 217   fHistosTrComb(0),
 218   fHistosTrCombU(0),
 219 //
 220   fHistosTrRcblPrim(0),
 221   fHistosTrRcblSec(0),
 222   fHistosCustom(0),
 223 //
 224   fEtaMin(-3.0),
 225   fEtaMax(3.0),
 226   fZVertexMin(-20),
 227   fZVertexMax( 20),
 228 //
 229   fScaleDTBySin2T(kFALSE),
 230   fCutOnDThetaX(kFALSE),
 231   fNStdDev(1.),
 232   fDPhiWindow(0.08),
 233   fDThetaWindow(0.025),
 234   fDPhiShift(0.0045),
 235   fPhiOverlapCut(0.005),
 236   fZetaOverlap(0.05),
 237   fPhiRot(0.),
 238   fInjScale(1.),
 239   fRemoveOverlaps(kFALSE),
 240 //
 241   fDPhiSCut(0.06),
 242   fNStdCut(1.),
 243   fMCV0Scale(0.7520),
 244 //
 245   fMultReco(0),
 246   fRPTree(0),
 247   fRPTreeMix(0),
 248   fStack(0),
 249   fMCEvent(0),
 250   //
 251   fNPart(0),
 252   fNBColl(0),
 253   fCurrCentBin(-1),
 254   fNCentBins(0),
 255   fUseCentralityVar(kCentV0M)
 256 {
 257   // Constructor
 258
 259   DefineOutput(1, TList::Class());
 260   //
 261   SetScaleDThetaBySin2T();
 262   SetNStdDev();
 263   SetPhiWindow();
 264   SetThetaWindow();
 265   SetPhiShift();
 266   SetPhiOverlapCut();
 267   SetZetaOverlapCut();
 268   SetPhiRot();
 269   SetRemoveOverlaps();
 270   //
 271 }
 272
 273 //________________________________________________________________________
 274 AliTrackletTaskMulti::~AliTrackletTaskMulti()
 275 {
 276   // Destructor
 277   // histograms are in the output list and deleted when the output
 278   // list is deleted by the TSelector dtor
 279   if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {  //RRR
 280     printf("Deleteing output\n");
 281     delete fOutput;
 282     fOutput = 0;
 283   }
 284   //
 285   delete fMultReco;
 286   //
 287   delete fHistosTrData;
 288   delete fHistosTrPrim;
 289   delete fHistosTrSec;
 290   delete fHistosTrComb;
 291   delete fHistosTrCombU;
 292   delete fHistosTrInj;
 293   delete fHistosTrRot;
 294   delete fHistosTrMix;
 295   delete fHistosTrRcblPrim;
 296   delete fHistosTrRcblSec;
 297   delete fHistosCustom;
 298   //
 299 }
 300
 301 //________________________________________________________________________
 302 void AliTrackletTaskMulti::UserCreateOutputObjects()
 303 {
 304   //
 305   OpenFile(1);
 306   fOutput = new TList();
 307   fOutput->SetOwner();
 308   //
 309   //
 310   Bool_t needGeom = GetDoNormalReco() || GetDoInjection() || GetDoRotation() || GetDoMixing();
 311   if (needGeom) {
 312     AliCDBManager *man = AliCDBManager::Instance();
 313     if (fUseMC) {
 314       Bool_t newGeom = kTRUE;
 315       man->SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Residual");
 316       if (newGeom) {
 317         // new geom
 318         AliCDBEntry*  obj = man->Get("GRP/Geometry/Data",130844,8);
 319         AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
 320         if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",130844,6,-1)) AliFatal("Failed to misalign geometry");
 321       }
 322       else {
 323         // old geom
 324         AliCDBEntry*  obj = man->Get("GRP/Geometry/Data",130845,7);
 325         AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
 326         if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",130845,5,-1)) AliFatal("Failed to misalign geometry");
 327       }
 328     }
 329     else {
 330       man->SetDefaultStorage("alien://Folder=/alice/data/2010/OCDB"); //man->SetRun(137366);
 331       AliCDBEntry*  obj = man->Get("GRP/Geometry/Data",137366);
 332       AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
 333       if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",137366,-1,-1)) AliFatal("Failed to misalign geometry");
 334     }
 335   }
 336     //
 337   // Create histograms
 338   fNCentBins = sizeof(fgkCentPerc)/sizeof(Float_t)-1;
 339   //---------------------------------------------Standard histos per tracklet type--->>
 340   UInt_t hPattern = 0xffffffff;
 341   fHistosTrData                      = BookHistosSet("TrData",hPattern);
 342   hPattern &= ~(BIT(kHEtaZvSPD1));  // fill single clusters for "data" only
 343   if (GetDoInjection()) fHistosTrInj = BookHistosSet("TrInj",hPattern);
 344   if (GetDoRotation())  fHistosTrRot = BookHistosSet("TrRot",hPattern);
 345   if (GetDoMixing())    fHistosTrMix = BookHistosSet("TrMix",hPattern);
 346   if (fUseMC) {
 347     fHistosTrPrim  = BookHistosSet("TrPrim",hPattern);
 348     fHistosTrSec   = BookHistosSet("TrSec",hPattern);
 349     fHistosTrComb  = BookHistosSet("TrComb",hPattern);
 350     fHistosTrCombU = BookHistosSet("TrCombU",hPattern);
 351     if (fCheckReconstructables) {
 352       fHistosTrRcblPrim = BookHistosSet("TrRcblPrim",hPattern);
 353       fHistosTrRcblSec  = BookHistosSet("TrRcblSec",hPattern);
 354     }
 355   }
 356   //---------------------------------------------Standard histos per tracklet type---<<
 357   //
 358   //---------------------------------------------Custom Histos----------------------->>
 359   // put here any non standard histos
 360   fHistosCustom = BookCustomHistos();
 361   //
 362   //---------------------------------------------Custom Histos-----------------------<<
 363   int nhist = fOutput->GetEntries();
 364   for (int i=0;i<nhist;i++) {
 365     TObject* hst = fOutput->At(i);
 366     if (!hst || !(hst->InheritsFrom(TH1::Class()))) continue;
 367     ((TH1*)hst)->Sumw2();
 368   }
 369   //
 370   if (fUseCentralityVar<0||fUseCentralityVar>kNCentTypes) {
 371     printf("Wrong centrality type %d\n",fUseCentralityVar);
 372     exit(1);
 373   }
 374   AliInfo(Form("Centrality type selected: %s\n",fgCentSelName[fUseCentralityVar]));
 375   PostData(1, fOutput);
 376   //
 377 }
 378
 379 //________________________________________________________________________
 380 void AliTrackletTaskMulti::UserExec(Option_t *)
 381 {
 382   // Main loop
 383   //
 384   static Bool_t statOK = kFALSE;
 385   if (!statOK) {
 386     RegisterStat();
 387     statOK = kTRUE;
 388   }
 389   //
 390   Bool_t needRecPoints = GetDoNormalReco() || GetDoInjection() || GetDoRotation() || GetDoMixing();
 391   //
 392   AliAnalysisManager* anMan = AliAnalysisManager::GetAnalysisManager();
 393   fRPTree = fRPTreeMix = 0;
 394   AliESDInputHandler *handler = (AliESDInputHandler*)anMan->GetInputEventHandler();
 395   AliESDInputHandlerRP *handRP = 0;
 396   if (needRecPoints) {
 397     handRP = (AliESDInputHandlerRP*)handler;
 398     if (!handRP) { printf("No RP handler\n"); return; }
 399   }
 400   AliESDEvent *esd  = handler->GetEvent();
 401   if (!esd) { printf("No AliESDEvent\n"); return; }
 402   //
 403   // do we need to initialize the field?
 404   AliMagF* field = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
 405   if (!field && !esd->InitMagneticField()) {printf("Failed to initialize the B field\n");return;}
 406   //
 407   //
 408   TH1* hstat = (TH1*)fHistosCustom->UncheckedAt(kHStat);
 409   hstat->Fill(kEvTot0); // RS
 410   const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
 411   const AliMultiplicity* multESD = esd->GetMultiplicity();
 412   //
 413   if (vtxESD->GetNContributors()<1) return;
 414   TString vtxTyp = vtxESD->GetTitle();
 415   if (vtxTyp.Contains("vertexer: Z")) {
 416     if (vtxESD->GetDispersion()>0.04) return;
 417     if (vtxESD->GetZRes()>0.25) return;
 418   }
 419   //
 420   AliCentrality *centrality = esd->GetCentrality();
 421   //
 422   hstat->Fill(kEvTot); // RS
 423   //
 424   Double_t esdvtx[3];
 425   vtxESD->GetXYZ(esdvtx);
 426   for (int i=3;i--;) fESDVtx[i] = esdvtx[i];
 427   //
 428   float vtxf[3] = {vtxESD->GetX(),vtxESD->GetY(),vtxESD->GetZ()};
 429   //
 430   // registed Ntracklets and ZVertex of the event
 431   ((TH1*)fHistosCustom->UncheckedAt(kHZVtxNoSel))->Fill(esdvtx[2]);
 432   //
 433   if(vtxf[2] < fZVertexMin || vtxf[2] > fZVertexMax) return;
 434   //
 435   //  centrality->Dump();
 436   Float_t multV0=0;
 437   AliESDVZERO* esdV0 = esd->GetVZEROData();
 438   if (esdV0) {
 439     multV0 = esdV0->GetMTotV0A()+esdV0->GetMTotV0C();
 440     if (fUseMC) multV0 *= fMCV0Scale;
 441   }
 442   ((TH1*)fHistosCustom->UncheckedAt(kHV0NoSel))->Fill(multV0);
 443   //
 444   float nSPD2 = multESD->GetNumberOfITSClusters(1);
 445   ((TH1*)fHistosCustom->UncheckedAt(kHNClSPD2NoSel))->Fill(nSPD2);
 446   //
 447   //------------------------------------------------------
 448   AliESDZDC *esdZDC = esd->GetESDZDC();
 449   float zdcEnergy=0,zemEnergy=0;
 450   if (esdZDC) {
 451     zdcEnergy = (esdZDC->GetZDCN1Energy() + esdZDC->GetZDCP1Energy() + esdZDC->GetZDCN2Energy()+ esdZDC->GetZDCP2Energy());
 452     zemEnergy = (esdZDC->GetZDCEMEnergy(0)+esdZDC->GetZDCEMEnergy(1))/8.;
 453   }
 454   ((TH2*)fHistosCustom->UncheckedAt(kHZDCZEMNoSel))->Fill(zemEnergy,zdcEnergy);
 455   //
 456   Float_t centPercentile = centrality->GetCentralityPercentileUnchecked(fgCentSelName[fUseCentralityVar]);
 457
 458   // temporary >>>>>>>>>>>>>>>>>>>>>>>>
 459   /*
 460   if (fUseCentralityVar==kCentZEMvsZDC) {
 461     float zdcEn = zdcEnergy;
 462     float zemEn = zemEnergy;
 463     Float_t slope;
 464     Float_t zdcPercentile;
 465     if (zemEn > 295.) {
 466       slope = (zdcEn + 15000.)/(zemEn - 295.);
 467       slope += 2.23660e+02;
 468       zdcPercentile = (TMath::ATan(slope) - 1.56664)/8.99571e-05;
 469       if (zdcPercentile<0) zdcPercentile = 0;
 470     }
 471     else zdcPercentile = 100;
 472     centPercentile = zdcPercentile;
 473   }
 474   else if (fUseMC && fUseCentralityVar==kCentV0M) {
 475     static TH1* fHtempV0M = 0;
 476     if (!fHtempV0M) {
 477       TString flcp = Form("%s/COMMON/CENTRALITY/data/pass%d/AliCentralityBy1D_%d.root", AliAnalysisManager::GetOADBPath(), 2, 139172);
 478       gSystem->ExpandPathName(flcp);
 479       TFile* flc  = TFile::Open(flcp);
 480       fHtempV0M  = (TH1*) (flc->Get("hmultV0_percentile"));
 481       fHtempV0M->SetDirectory(0);
 482       flc->Close();
 483       delete flc;
 484       printf("Getting V0 centrality histo from %s -> %p\n",flcp.Data(),fHtempV0M);
 485     }
 486     if (!fHtempV0M) AliFatal("No V0 centrality histo");
 487     centPercentile = fHtempV0M->GetBinContent(fHtempV0M->FindBin((multV0)));
 488   }
 489   */
 490   // temporary >>>>>>>>>>>>>>>>>>>>>>>>
 491
 492   fCurrCentBin = GetCentralityBin(centPercentile);
 493   //
 494   printf("CentPerc: %f : Bin %d ZDC: %f ZEM: %f\n",centPercentile, fCurrCentBin, zdcEnergy,zemEnergy);
 495   if (fCurrCentBin<0) {
 496     //printf("Reject: %.1f : V0:%.1f V0Cor:%.1f V0CR:%.1f SPD2c:%.1f\n",mltTst, multV0,multV0Corr,multV0CorrResc,nSPD2Corr);
 497     return;
 498   }
 499   //
 500   ((TH1*)fHistosCustom->UncheckedAt(kHStatCentBin))->Fill(fCurrCentBin);
 501   ((TH1*)fHistosCustom->UncheckedAt(kHStatCent))->Fill(centPercentile);
 502   //
 503   AliMCEventHandler* eventHandler = 0;
 504   fMCEvent = 0;
 505   fStack = 0;
 506   //
 507   if (fUseMC) {
 508     eventHandler = (AliMCEventHandler*)anMan->GetMCtruthEventHandler();
 509     if (!eventHandler) { printf("ERROR: Could not retrieve MC event handler\n"); return; }
 510     fMCEvent = eventHandler->MCEvent();
 511     if (!fMCEvent) { printf("ERROR: Could not retrieve MC event\n"); return; }
 512     fStack = fMCEvent->Stack();
 513     if (!fStack) { printf("Stack not available\n"); return; }
 514   }
 515   //
 516   if (needRecPoints) {
 517     fRPTree = handRP->GetTreeR("ITS");
 518     if (!fRPTree) { AliError(" Invalid ITS cluster tree !\n"); return; }
 519   }
 520   //
 521   // =============================================================================>>>
 522   // MC Generator info
 523   AliGenEventHeader* mcGenH = 0;
 524   fNPart  = 0;
 525   fNBColl = 0;
 526   if (fUseMC) {
 527     mcGenH = fMCEvent->GenEventHeader();
 528     if (mcGenH->InheritsFrom(AliGenHijingEventHeader::Class())) {
 529       AliGenHijingEventHeader* hHijing = (AliGenHijingEventHeader*)mcGenH;
 530       fNPart  = (hHijing->ProjectileParticipants()+hHijing->TargetParticipants())/2.;
 531       fNBColl = hHijing->NN()+hHijing->NNw()+hHijing->NwN()+hHijing->NwNw();
 532     }
 533     else if (mcGenH->InheritsFrom(AliGenDPMjetEventHeader::Class())) {
 534       AliGenDPMjetEventHeader* hDpmJet = (AliGenDPMjetEventHeader*)mcGenH;
 535       fNPart  = (hDpmJet->ProjectileParticipants()+hDpmJet->TargetParticipants())/2.;
 536       fNBColl = hDpmJet->NN()+hDpmJet->NNw()+hDpmJet->NwN()+hDpmJet->NwNw();
 537     }
 538     else {} // unknown generator
 539   }
 540   //
 541   // register Ntracklets and ZVertex of the event
 542   ((TH2*)fHistosCustom->UncheckedAt(kHZVtx))->Fill(esdvtx[2],fCurrCentBin);
 543   ((TH2*)fHistosCustom->UncheckedAt(kHV0))->Fill(multV0,fCurrCentBin);
 544   ((TH2*)fHistosCustom->UncheckedAt(kHNClSPD2))->Fill(nSPD2,fCurrCentBin);
 545   ((TH3*)fHistosCustom->UncheckedAt(kHZDCZEM))->Fill(zemEnergy,zdcEnergy,fCurrCentBin);
 546   //
 547   if (fUseMC) FillMCPrimaries();
 548   //
 549   // normal reconstruction
 550   hstat->Fill(kBinEntries+kEvProcData + kEntriesPerBin*fCurrCentBin);
 551   //
 552   if (GetDoNormalReco() || GetDoInjection()) { // for the injection the normal reco should be done
 553     InitMultReco();
 554     fMultReco->Run(fRPTree, vtxf);
 555     printf("Multiplicity Reconstructed:\n");
 556     AliMultiplicity* mlt = fMultReco->GetMultiplicity();
 557     if (mlt) mlt->Print();
 558     if (GetDoNormalReco()) FillHistos(kData,mlt);
 559     FillClusterInfo();
 560     //
 561   }
 562   if (!GetDoNormalReco()) {
 563     FillHistos(kData,multESD); // fill data histos from ESD
 564     FillClusterInfoFromMult(multESD, vtxf[2] );
 565   }
 566   //
 567   // Injection: it must come right after the normal reco since needs its results
 568   if (GetDoInjection()) {
 569     if (!fMultReco) InitMultReco(); // in principle, not needed, the reco is created above
 570     fMultReco->SetRecType(AliITSMultRecBg::kBgInj);
 571     fMultReco->Run(fRPTree, vtxf);
 572     printf("Multiplicity from Injection:\n");
 573     AliMultiplicity* mlt = fMultReco->GetMultiplicity();
 574     if (mlt) mlt->Print();
 575     hstat->Fill(kBinEntries + kEvProcInj + kEntriesPerBin*fCurrCentBin);
 576     FillHistos(kBgInj,mlt);
 577   }
 578   //
 579   // Rotation
 580   if (GetDoRotation()) {
 581     InitMultReco();
 582     fMultReco->SetRecType(AliITSMultRecBg::kBgRot);
 583     fMultReco->SetPhiRotationAngle(fPhiRot);
 584     fMultReco->Run(fRPTree, vtxf);
 585     printf("Multiplicity from Rotation:\n");
 586     AliMultiplicity* mlt = fMultReco->GetMultiplicity();
 587     if (mlt) mlt->Print();
 588     hstat->Fill(kBinEntries + kEvProcRot + kEntriesPerBin*fCurrCentBin);
 589     FillHistos(kBgRot,mlt);
 590   }
 591   //
 592   if (GetDoMixing()) {
 593     /*
 594     AliMixEventInputHandler* handToMix = (AliMixEventInputHandler*)handRP->MixingHandler();
 595     if (!handToMix) { printf("No Mixing handler\n"); return; }
 596     handToMix->GetEntry();
 597     if(handToMix->MixedEventNumber()<1) {printf("Mixing: No enough events in pool\n"); return;}
 598     AliESDInputHandlerRP* handRPMix = (AliESDInputHandlerRP*) handToMix->InputEventHandler(0);
 599
 600     if (!handRPMix) { printf("No Mixing RP handler\n"); return; }
 601     fRPTreeMix = handRPMix->GetTreeR("ITS");
 602     if (!fRPTreeMix) { AliError(" Invalid ITS cluster tree of the 2nd event!\n"); return; }
 603     //
 604     AliESDEvent *esdMix = handRPMix->GetEvent();
 605     const AliESDVertex* vtxESDmix = esdMix->GetVertex();
 606     ((TH2*)fHistosCustom->UncheckedAt(kHZVtxMixDiff))->Fill(vtxESDmix->GetZ()-esdvtx[2],fCurrCentBin);
 607     ((TH2*)fHistosCustom->UncheckedAt(kHNTrMixDiff) )->
 608       Fill(esdMix->GetMultiplicity()->GetNumberOfTracklets() - multESD->GetNumberOfTracklets(),fCurrCentBin);
 609     //
 610     InitMultReco();
 611     fMultReco->SetRecType(AliITSMultRecBg::kBgMix);
 612     fMultReco->Run(fRPTree, vtxf,fRPTreeMix);
 613     printf("Multiplicity from Mixing:\n");
 614     AliMultiplicity* mlt = fMultReco->GetMultiplicity();
 615     if (mlt) mlt->Print();
 616     hstat->Fill(kBinEntries + kEvProcMix + kEntriesPerBin*fCurrCentBin);
 617     FillHistos(kBgMix,mlt);
 618     */
 619     //
 620   }
 621   // =============================================================================<<<
 622   //
 623   if (fMultReco) delete fMultReco;
 624   fMultReco = 0;
 625   //
 626 }
 627
 628
 629 //________________________________________________________________________
 630 void AliTrackletTaskMulti::RegisterStat()
 631 {
 632   TH1* hstat;
 633   TList *lst = dynamic_cast<TList*>(GetOutputData(1));
 634   if (lst && (hstat=(TH1*)lst->FindObject("hStat"))) {
 635     Info("Terminate","Registering used settings");
 636     // fill used settings
 637     hstat->Fill(kOneUnit,1.);
 638     hstat->Fill(kCentVar,fUseCentralityVar);
 639     hstat->Fill(kDPhi,fDPhiWindow);
 640     hstat->Fill(kDTht,fDThetaWindow);
 641     hstat->Fill(kNStd,fNStdDev);
 642     hstat->Fill(kPhiShift,fDPhiShift);
 643     hstat->Fill(kThtS2,fScaleDTBySin2T);
 644     hstat->Fill(kThtCW,fCutOnDThetaX);
 645     hstat->Fill(kPhiOvl,fPhiOverlapCut);
 646     hstat->Fill(kZEtaOvl,fZetaOverlap);
 647     hstat->Fill(kNoOvl,fRemoveOverlaps);
 648     //
 649     hstat->Fill(kPhiRot,fPhiRot);
 650     hstat->Fill(kInjScl,fInjScale);
 651     hstat->Fill(kEtaMin,fEtaMin);
 652     hstat->Fill(kEtaMax,fEtaMax);
 653     hstat->Fill(kZVMin,fZVertexMin);
 654     hstat->Fill(kZVMax,fZVertexMax);
 655     //
 656     hstat->Fill(kDPiSCut,fDPhiSCut);
 657     hstat->Fill(kNStdCut,fNStdCut);
 658     hstat->Fill(kMCV0Scale, fMCV0Scale);
 659     //
 660   }
 661   //
 662 }
 663
 664
 665 //________________________________________________________________________
 666 void AliTrackletTaskMulti::Terminate(Option_t *)
 667 {
 668   Printf("Terminating...");
 669   RegisterStat();
 670 }
 671
 672
 673 //_________________________________________________________________________
 674 void AliTrackletTaskMulti::InitMultReco()
 675 {
 676   // create mult reconstructor
 677   if (fMultReco) delete fMultReco;
 678   fMultReco = new AliITSMultRecBg();
 679   fMultReco->SetCreateClustersCopy(kTRUE);
 680   fMultReco->SetScaleDThetaBySin2T(fScaleDTBySin2T);
 681   fMultReco->SetNStdDev(fNStdDev);
 682   fMultReco->SetPhiWindow( fDPhiWindow );
 683   fMultReco->SetThetaWindow( fDThetaWindow );
 684   fMultReco->SetPhiShift( fDPhiShift );
 685   fMultReco->SetRemoveClustersFromOverlaps(fRemoveOverlaps);
 686   fMultReco->SetPhiOverlapCut(fPhiOverlapCut);
 687   fMultReco->SetZetaOverlapCut(fZetaOverlap);
 688   fMultReco->SetHistOn(kFALSE);
 689   fMultReco->SetRecType( AliITSMultRecBg::kData );
 690 }
 691
 692 //_________________________________________________________________________
 693 TObjArray* AliTrackletTaskMulti::BookCustomHistos()
 694 {
 695   // book custom histos, not related to specific tracklet type
 696   TObjArray* histos = new TObjArray();
 697   TH1F* hstat;
 698   //
 699   // ------------ job parameters, statistics ------------------------------>>>
 700   int nbs = kBinEntries + fNCentBins*kEntriesPerBin;
 701   hstat = new TH1F("hStat","Run statistics",nbs,0.5,nbs+0.5);
 702   //
 703   hstat->GetXaxis()->SetBinLabel(kEvTot0, "Ev.Tot0");
 704   hstat->GetXaxis()->SetBinLabel(kEvTot, "Ev.Tot");
 705   hstat->GetXaxis()->SetBinLabel(kOneUnit,"ScaleMerge");
 706   hstat->GetXaxis()->SetBinLabel(kNWorkers,"Workers");
 707   //
 708   hstat->GetXaxis()->SetBinLabel(kCentVar,  fgCentSelName[fUseCentralityVar]);
 709   hstat->GetXaxis()->SetBinLabel(kDPhi,  "#Delta#varphi");
 710   hstat->GetXaxis()->SetBinLabel(kDTht,  "#Delta#theta");
 711   hstat->GetXaxis()->SetBinLabel(kNStd,  "N.std");
 712   hstat->GetXaxis()->SetBinLabel(kPhiShift,"#delta#varphi");
 713   hstat->GetXaxis()->SetBinLabel(kThtS2,"scale #Delta#theta");
 714   hstat->GetXaxis()->SetBinLabel(kPhiOvl,"#varpho_{Ovl}");
 715   hstat->GetXaxis()->SetBinLabel(kZEtaOvl,"#z_{Ovl}");
 716   hstat->GetXaxis()->SetBinLabel(kNoOvl, "rem.ovl");
 717   //
 718   hstat->GetXaxis()->SetBinLabel(kPhiRot,"#varphi_{rot}");
 719   hstat->GetXaxis()->SetBinLabel(kInjScl,"inj");
 720   hstat->GetXaxis()->SetBinLabel(kEtaMin,"#eta_{min}");
 721   hstat->GetXaxis()->SetBinLabel(kEtaMax,"#eta_{max}");
 722   hstat->GetXaxis()->SetBinLabel(kZVMin,"ZV_{min} cut");
 723   hstat->GetXaxis()->SetBinLabel(kZVMax,"ZV_{max} cut");
 724   //
 725   hstat->GetXaxis()->SetBinLabel(kDPiSCut,"#Delta#varphi-#delta_{#phi} cut");
 726   hstat->GetXaxis()->SetBinLabel(kNStdCut,"#Delta cut");
 727   //
 728   hstat->GetXaxis()->SetBinLabel(kMCV0Scale,"MC V0 scale");
 729   //
 730   for (int i=0;i<fNCentBins;i++) {
 731     TString bnt = "b"; bnt+= i;
 732     int offs = kBinEntries + i*kEntriesPerBin;
 733     hstat->GetXaxis()->SetBinLabel(offs + kEvProcData, bnt+" Ev.ProcData");
 734     hstat->GetXaxis()->SetBinLabel(offs + kEvProcInj,  bnt+" Ev.ProcInj");
 735     hstat->GetXaxis()->SetBinLabel(offs + kEvProcRot,  bnt+" Ev.ProcRot");
 736     hstat->GetXaxis()->SetBinLabel(offs + kEvProcMix,  bnt+" Ev.ProcMix");
 737     //
 738   }
 739   //
 740   hstat->Fill(kNWorkers);
 741   //
 742   AddHisto(histos,hstat,kHStat);
 743   //
 744   // ------------------------ events per centrality bin ----------------------
 745   TH1D* hCentAx = new TH1D("EvCentr","Events per centrality",fNCentBins,fgkCentPerc);
 746   hCentAx->GetXaxis()->SetTitle("Centrality parameter");
 747   AddHisto(histos,hCentAx,kHStatCent);
 748   //
 749   TH1D* hCentBin = new TH1D("EvCentrBin","Events per centrality bin",fNCentBins,-0.5,fNCentBins-0.5);
 750   hCentBin->GetXaxis()->SetTitle("Centrality Bin");
 751   AddHisto(histos,hCentBin,kHStatCentBin);
 752   //
 753   // ------------ job parameters, statistics ------------------------------<<<
 754   //
 755   //  double etaMn=-3,etaMx=3;
 756   double zMn=-30, zMx=30;
 757   //  int nEtaBins = int((etaMx-etaMn)/kEtaBinWidth);
 758   //  if (nEtaBins<1) nEtaBins = 1;
 759   //
 760   int nZVBins = int(zMx-zMn);
 761   if (nZVBins<1) nZVBins = 1;
 762   //
 763   // Z vertex distribution for events before selection
 764   TH1F* hzvns = new  TH1F("zvNoSel","Z vertex before selection",nZVBins,zMn,zMx);
 765   hzvns->GetXaxis()->SetTitle("Zvertex");
 766   AddHisto(histos,hzvns,kHZVtxNoSel);
 767   //
 768   // V0 for events before selection
 769   int nbmltV0 = 250;
 770   double maxmltV0 = 25000;
 771   //
 772   TH1F* hnV0ns = new  TH1F("V0NoSel","V0 signal Before Cent. Selection",nbmltV0,0,maxmltV0);
 773   hnV0ns->GetXaxis()->SetTitle("V0 signal");
 774   AddHisto(histos,hnV0ns,kHV0NoSel);
 775   //
 776   // N SPD2 clusters
 777   int nbmltSPD2 = 175;
 778   double maxmltSPD2 = 7000;
 779   TH1F* hncl2ns = new  TH1F("NClustersSPD2NoSel","N Clusters on SPD2 Before Cent Selection",nbmltSPD2,0,maxmltSPD2);
 780   hncl2ns->GetXaxis()->SetTitle("N Clus SPD2");
 781   AddHisto(histos,hncl2ns,kHNClSPD2NoSel);
 782   //
 783   int nbzdc=50;
 784   double maxZDC=6000., maxZEM=2500.;
 785   TH2F* hzdczemns = new  TH2F("ZDCZEMNoSel","ZDC vs ZEM Before Cent Selection",
 786                               nbzdc,0,maxZEM,nbzdc,0,maxZDC);
 787   hzdczemns->GetXaxis()->SetTitle("ZEM");
 788   hzdczemns->GetXaxis()->SetTitle("ZDC");
 789   AddHisto(histos,hzdczemns,kHZDCZEMNoSel);
 790   //
 791   TH2F* hzv = new  TH2F("zv","Z vertex after Selection per Cent.Bin",nZVBins,zMn,zMx, fNCentBins, -0.5,fNCentBins-0.5);
 792   hzv->GetXaxis()->SetTitle("Zvertex");
 793   hzv->GetYaxis()->SetTitle("Cent.Bin ID");
 794   AddHisto(histos,hzv,kHZVtx);
 795   //
 796   // V0
 797   TH2F* hnV0 = new  TH2F("V0","V0 signal per Cent.Bin ",nbmltV0,0,maxmltV0, fNCentBins, -0.5,fNCentBins-0.5);
 798   hnV0->GetXaxis()->SetTitle("V0 signal");
 799   hnV0->GetYaxis()->SetTitle("Cent.Bin ID");
 800   AddHisto(histos,hnV0,kHV0);
 801   //
 802   // N SPD2 clusters
 803   TH2F* hncl2 = new  TH2F("NClustersSPD2","N Clusters on SPD2 per Cent.Bin",nbmltSPD2,0,maxmltSPD2, fNCentBins, -0.5,fNCentBins-0.5);
 804   hncl2->GetXaxis()->SetTitle("N Clus SPD2");
 805   hncl2->GetYaxis()->SetTitle("Cent.Bin ID");
 806   AddHisto(histos,hncl2,kHNClSPD2);
 807   //
 808   // ZDCZEM
 809   TH3F* hzdczem = new TH3F("ZDCZEM","ZDC vs ZEM per Cent.Bin",nbzdc,0,maxZEM,nbzdc,0,maxZDC, fNCentBins, -0.5,fNCentBins-0.5);
 810   hzdczem->GetXaxis()->SetTitle("ZEM");
 811   hzdczem->GetYaxis()->SetTitle("ZDC");
 812   AddHisto(histos,hzdczem,kHZDCZEM);
 813   //
 814   //----------------------------------------------------------------------
 815   int nEtaBinsS = int((fEtaMax-fEtaMin)/kEtaBinWidth);
 816   if (nEtaBinsS<1) nEtaBinsS = 1;
 817   //
 818   int nZVBinsS = int(fZVertexMax-fZVertexMin);
 819   if (nZVBinsS<1) nZVBinsS = 1;
 820
 821   if (fUseMC) {
 822     // Z vertex vs Eta distribution for primaries
 823     char buffn[100],bufft[500];
 824     for (int ib=0;ib<fNCentBins;ib++) {
 825       sprintf(buffn,"b%d_zvEtaPrimMC",ib);
 826       sprintf(bufft,"bin%d Zvertex vs #eta PrimMC",ib);
 827       TH2F* hzvetap = new  TH2F(buffn,bufft, nEtaBinsS,fEtaMin,fEtaMax,nZVBinsS,fZVertexMin,fZVertexMax);
 828       hzvetap->GetXaxis()->SetTitle("#eta");
 829       hzvetap->GetYaxis()->SetTitle("Zvertex");
 830       AddHisto(histos,hzvetap,kHZVEtaPrimMC+ib);
 831     }
 832     //
 833     // <n> primaries according to MC generator
 834     TH1F* hnprimM = new  TH1F("nPrimMean","<N> primaries",fNCentBins, -0.5,fNCentBins-0.5);
 835     hnprimM->GetXaxis()->SetTitle("Cent.Bin ID");
 836     AddHisto(histos,hnprimM,kHNPrimMeanMC);
 837     //
 838     // <n> primaries per part.pair according to MC generator
 839     TH1F* hnprim2partM = new  TH1F("nPrim2Part","<N> primaries per part.pair",fNCentBins, -0.5,fNCentBins-0.5);
 840     hnprim2partM->GetXaxis()->SetTitle("Cent.Bin ID");
 841     AddHisto(histos,hnprim2partM,kHNPrim2PartMC);
 842     //
 843     // <n> primaries per part.pair vs npart.pair according to MC generator
 844     TH2F* hnprim2partNp = new  TH2F("nPrim2Part_vs_NPart","<N> primaries per part.pair vs N part.pairs",105,0,210,200,0,40);
 845     hnprim2partNp->GetXaxis()->SetTitle("N.part.pairs");
 846     hnprim2partNp->GetYaxis()->SetTitle("N.prim/N.part.pairs");
 847     AddHisto(histos,hnprim2partNp,kHNPrim2PartNpMC);
 848     //
 849     // <n> primaries per b.coll vs npart.pair according to MC generator
 850     TH2F* hnprim2BCollNp = new  TH2F("nPrim2BColl_vs_NPart","<N> primaries per bin.coll vs N part.pairs",105,0,210,200,0,40);
 851     hnprim2BCollNp->GetXaxis()->SetTitle("N.part.pairs");
 852     hnprim2BCollNp->GetYaxis()->SetTitle("N.prim/N.bin.coll.");
 853     AddHisto(histos,hnprim2BCollNp,kHNPrim2BCollNpMC);
 854     //
 855     // <n> primaries per bin.coll. according to MC generator
 856     TH1F* hnprim2BCollM = new  TH1F("nPrim2BColl","<N> primaries per bin.coll",fNCentBins, -0.5,fNCentBins-0.5);
 857     hnprim2BCollM->GetXaxis()->SetTitle("Cent.Bin ID");
 858     AddHisto(histos,hnprim2BCollM,kHNPrim2BCollMC);
 859     //
 860     // n participants according to MC generator
 861     TH2F* hnpart = new  TH2F("nPart","N participant pairs",210,0,210,fNCentBins, -0.5,fNCentBins-0.5);
 862     hnpart->GetXaxis()->SetTitle("N part. pairs");
 863     hnpart->GetYaxis()->SetTitle("Cent.Bin ID");
 864     AddHisto(histos,hnpart,kHNPartMC);
 865     //
 866     // <n> participants according to MC generator
 867     TH1F* hnpartM = new  TH1F("nPartMean","<N> participant pairs",fNCentBins, -0.5,fNCentBins-0.5);
 868     hnpartM->GetXaxis()->SetTitle("Cent.Bin ID");
 869     AddHisto(histos,hnpartM,kHNPartMeanMC);
 870     //
 871     // n bin coll. according to MC generator
 872     TH2F* hnbcoll = new  TH2F("nBColl","N bin. coll",2000,0,2000,fNCentBins, -0.5,fNCentBins-0.5);
 873     hnbcoll->GetXaxis()->SetTitle("N bin. coll");
 874     hnbcoll->GetYaxis()->SetTitle("Cent.Bin ID");
 875     AddHisto(histos,hnbcoll,kHNBCollMC);
 876     //
 877     // <n> bin col according to MC generator
 878     TH1F* hnbcollM = new  TH1F("nBCollMean","<N> bin.colls",fNCentBins, -0.5,fNCentBins-0.5);
 879     hnbcollM->GetXaxis()->SetTitle("Cent.Bin ID");
 880     AddHisto(histos,hnbcollM,kHNBCollMeanMC);
 881     //
 882   }
 883   //
 884   if (GetDoMixing()) {
 885     //
 886     // Difference in Z vertex for mixed events
 887     TH2F* hzdiff = new TH2F("MixSPDVertexDiff","SPD #Delta Z Vertex distribution per mult bin ",100,-5,5, fNCentBins, -0.5,fNCentBins-0.5);
 888     hzdiff->GetXaxis()->SetTitle("#Delta Z Vertex [cm]");
 889     hzdiff->GetYaxis()->SetTitle(Form("Entries / %1.2f [cm] per mult bin",10./100.));
 890     AddHisto(histos,hzdiff,kHZVtxMixDiff);
 891     //
 892     // Difference in N tracklets for mixed events
 893     TH2F* hntdiff = new TH2F("MixNTrackletsDiff"," SPD tracklets Diff ",200,-1000,1000, fNCentBins, -0.5,fNCentBins-0.5);
 894     hntdiff->GetXaxis()->SetTitle("# tracklet diff");
 895     AddHisto(histos,hntdiff,kHNTrMixDiff);
 896   }
 897   //
 898   // --------------------------------------------------
 899   if (fUseMC) {
 900     int npdg = sizeof(fgkPDGNames)/sizeof(char*);
 901     TH2F* hpdgP = new TH2F("pdgPrim","primary PDG",npdg,0,npdg,fNCentBins, -0.5,fNCentBins-0.5);
 902     AddHisto(histos,hpdgP,kHPrimPDG);
 903     TH2F* hpdgS = new TH2F("pdgSec","secondary PDG",npdg,0,npdg,fNCentBins, -0.5,fNCentBins-0.5);
 904     AddHisto(histos,hpdgS,kHSecPDG);
 905     TH2F* hpdgPP = new TH2F("pdgPrimPar","primary parent PDG ",npdg,0,npdg,fNCentBins, -0.5,fNCentBins-0.5);
 906     AddHisto(histos,hpdgPP,kHPrimParPDG);
 907     TH2F* hpdgSP = new TH2F("pdgSecPar","secondary parent PDG",npdg,0,npdg,fNCentBins, -0.5,fNCentBins-0.5);
 908     AddHisto(histos,hpdgSP,kHSecParPDG);
 909     for (int i=0;i<npdg;i++) {
 910       hpdgP->GetXaxis()->SetBinLabel(i+1,fgkPDGNames[i]);
 911       hpdgS->GetXaxis()->SetBinLabel(i+1,fgkPDGNames[i]);
 912       hpdgPP->GetXaxis()->SetBinLabel(i+1,fgkPDGNames[i]);
 913       hpdgSP->GetXaxis()->SetBinLabel(i+1,fgkPDGNames[i]);
 914     }
 915   }
 916   //
 917   // -------------------------------------------------
 918   TH2F* hclinf=0;
 919   hclinf = new TH2F("cl0InfoUsed","#phi vs Z of used clusters, Lr0",64,-16,16, 80,0,2*TMath::Pi());
 920   AddHisto(histos,hclinf,kHClUsedInfoL0);
 921   hclinf = new TH2F("cl1InfoUsed","#phi vs Z of used clusters, Lr1",64,-16,16, 2*80,0,2*TMath::Pi());
 922   AddHisto(histos,hclinf,kHClUsedInfoL1);
 923   hclinf = new TH2F("cl0InfoAll","#phi vs Z of all clusters, Lr0",64,-16,16, 80,0,2*TMath::Pi());
 924   AddHisto(histos,hclinf,kHClAllInfoL0);
 925   hclinf = new TH2F("cl1InfoAll","#phi vs Z of all clusters, Lr1",64,-16,16, 2*80,0,2*TMath::Pi());
 926   AddHisto(histos,hclinf,kHClAllInfoL1);
 927   //
 928   // -------------------------------------------------
 929   histos->SetOwner(kFALSE);
 930   //
 931   return histos;
 932 }
 933
 934 //_________________________________________________________________________
 935 TObjArray* AliTrackletTaskMulti::BookHistosSet(const char* pref, UInt_t selHistos)
 936 {
 937   // book standard set of histos attaching the pref in front of the name/title
 938   //
 939   const int kNDPhiBins = 100;
 940   const int kNDThtBins = 100;
 941   int nDistBins = int(fNStdDev)*5;
 942   //
 943   int nEtaBins = int((fEtaMax-fEtaMin)/kEtaBinWidth);
 944   if (nEtaBins<1) nEtaBins = 1;
 945   //
 946   int nZVBins = int(fZVertexMax-fZVertexMin);
 947   if (nZVBins<1) nZVBins = 1;
 948   float dphir = fDPhiWindow*TMath::Sqrt(fNStdDev);
 949   float dthtr = fDThetaWindow*TMath::Sqrt(fNStdDev);
 950   //
 951   TObjArray* histos = new TObjArray();
 952   TH2F* h2;
 953   TH1F* h1;
 954   char buffn[100],bufft[500];
 955   //
 956   for (int ib=0;ib<fNCentBins;ib++) {
 957     //
 958     int offs = ib*kNStandardH;
 959     if (selHistos & (0x1<<kHEtaZvCut) ) {
 960       sprintf(buffn,"b%d_%s_ZvEtaCutT",ib,pref);
 961       sprintf(bufft,"bin%d (%s) Zv vs Eta with tracklet cut",ib,pref);
 962       h2 = new TH2F(buffn,bufft,nEtaBins,fEtaMin,fEtaMax, nZVBins, fZVertexMin,fZVertexMax);
 963       h2->GetXaxis()->SetTitle("#eta");
 964       h2->GetYaxis()->SetTitle("Zv");
 965       AddHisto(histos,h2,offs+kHEtaZvCut);
 966     }
 967     //
 968     if (selHistos & (0x1<<kHDPhiDTheta) ) {
 969       sprintf(buffn,"b%d_%s_dPhidTheta",ib,pref);
 970       sprintf(bufft,"bin%d (%s) #Delta#theta vs #Delta#varphi",ib,pref);
 971       h2 = new TH2F(buffn,bufft,kNDPhiBins,-dphir,dphir,kNDThtBins,-dthtr,dthtr);
 972       h2->GetXaxis()->SetTitle("#Delta#varphi [rad]");
 973       h2->GetYaxis()->SetTitle("#Delta#theta [rad]");
 974       AddHisto(histos,h2,offs+kHDPhiDTheta);
 975     }
 976     //
 977     if (selHistos & (0x1<<kHDPhiSDThetaX) ) {
 978       sprintf(buffn,"b%d_%s_dPhiSdThetaX",ib,pref);
 979       sprintf(bufft,"bin%d (%s) #Delta#theta%s vs #Delta#varphi-#delta_{#varphi}",ib,pref,fScaleDTBySin2T ? "/sin^{2}(#theta)":"");
 980       h2 = new TH2F(buffn,bufft,kNDPhiBins,-dphir,dphir,kNDThtBins,-dthtr,dthtr);
 981       h2->GetXaxis()->SetTitle("#Delta#varphi-#delta_{#varphi} [rad]");
 982       sprintf(bufft,"#Delta#theta%s",fScaleDTBySin2T ? "/sin^{2}(#theta)":"");
 983       h2->GetYaxis()->SetTitle(bufft);
 984       AddHisto(histos,h2,offs+kHDPhiSDThetaX);
 985     }
 986     //
 987     if (selHistos & (0x1<<kHWDist) ) {
 988       sprintf(buffn,"b%d_%s_WDist",ib,pref);
 989       sprintf(bufft,"bin%d #Delta=[(#Delta#varphi-#delta_{#varphi})/#sigma#varphi]^{2}+"
 990               "[#Delta#theta%s/#sigma#theta]^{2}",ib,fScaleDTBySin2T ? "*sin^{-2}(#theta)":"");
 991       h1 = new TH1F(buffn,bufft,nDistBins,0,fNStdDev);
 992       sprintf(bufft,"#Delta=[(#Delta#varphi-#delta_{#varphi})/#sigma#varphi]^{2}+"
 993               "[#Delta#theta%s/#sigma#theta]^{2}",fScaleDTBySin2T ? "*sin^{-2}(#theta)":"");
 994       h1->GetXaxis()->SetTitle(bufft);
 995       AddHisto(histos,h1,offs+kHWDist);
 996     }
 997     //
 998     if (selHistos & (0x1<<kHEtaZvSPD1) ) {
 999       sprintf(buffn,"b%d_%s_ZvEtaSPD1",ib,pref);
1000       sprintf(bufft,"bin%d (%s) Zv vs Eta SPD1 clusters",ib,pref);
1001       h2 = new TH2F(buffn,bufft,nEtaBins,fEtaMin,fEtaMax, nZVBins, fZVertexMin,fZVertexMax);
1002       h2->GetXaxis()->SetTitle("#eta");
1003       h2->GetYaxis()->SetTitle("Zv");
1004       AddHisto(histos,h2,offs+kHEtaZvSPD1);
1005     }
1006     //
1007   }
1008   //
1009   histos->SetOwner(kFALSE);
1010   return histos;
1011 }
1012
1013 //_________________________________________________________________________
1014 void AliTrackletTaskMulti::AddHisto(TObjArray* histos, TObject* h, Int_t at)
1015 {
1016   // add single histo to the set
1017   if (at>=0) histos->AddAtAndExpand(h,at);
1018   else       histos->Add(h);
1019   fOutput->Add(h);
1020 }
1021
1022 //_________________________________________________________________________
1023 void AliTrackletTaskMulti::FillHistos(Int_t type, const AliMultiplicity* mlt)
1024 {
1025   // fill histos of given type
1026   if (!mlt) return;
1027   //
1028   TObjArray* histos = 0;
1029   if      (type == kData)  histos = fHistosTrData;
1030   else if (type == kBgInj) histos = fHistosTrInj;
1031   else if (type == kBgRot) histos = fHistosTrRot;
1032   else if (type == kBgMix) histos = fHistosTrMix;
1033   //
1034   Bool_t fillMC = (type==kData) && fUseMC && fStack;
1035   //
1036   //
1037   //---------------------------------------- CHECK ------------------------------>>>
1038   TArrayF vtxMC;
1039   AliGenHijingEventHeader* pyHeader = 0;
1040   //
1041   if (fUseMC) {
1042     pyHeader = (AliGenHijingEventHeader*) fMCEvent->GenEventHeader();//header->GenEventHeader();
1043     pyHeader->PrimaryVertex(vtxMC);
1044   }
1045   //---------------------------------------- CHECK ------------------------------<<<
1046   //
1047   if (!histos) return;
1048   int ntr = mlt->GetNumberOfTracklets();
1049   for (int itr=ntr;itr--;) {
1050     //
1051     //---------------------------------------- CHECK ------------------------------>>>
1052     /*
1053     if (fUseMC) {
1054       Bool_t reject = kFALSE;
1055       while(1) {
1056         int lab0 = mlt->GetLabel(itr,0);
1057         int lab1 = mlt->GetLabel(itr,1);
1058         if (lab0!=lab1) break;
1059         if (!fStack->IsPhysicalPrimary(lab0)) break;
1060         //
1061         TParticle* part = fStack->Particle(lab0);
1062         Float_t dz = part->Vz() - vtxMC[2];
1063         if (TMath::Abs(dz)<1e-6) break;
1064         reject = kTRUE;
1065         break;
1066       }
1067       if (reject) continue;
1068     }
1069     */
1070     //---------------------------------------- CHECK ------------------------------<<<
1071     //
1072     double theta  = mlt->GetTheta(itr);
1073     double eta    = -TMath::Log(TMath::Tan(theta/2));
1074     if (eta<fEtaMin || eta>fEtaMax) continue;
1075     //
1076     double dtheta = mlt->GetDeltaTheta(itr);
1077     double dThetaX = dtheta;
1078     if (fScaleDTBySin2T) {
1079       double sint   =  TMath::Sin(theta);
1080       dThetaX /= (sint*sint);
1081     }
1082     if (fCutOnDThetaX && TMath::Abs(dThetaX)>fDThetaWindow) continue;
1083     //    double phi    = mlt->GetPhi(itr);
1084     double dphi   = mlt->GetDeltaPhi(itr);
1085     double dist   = mlt->CalcDist(itr);
1086     //
1087     FillHistosSet(histos,eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist);
1088     // special handling for mc info
1089     if (fillMC && fStack) {
1090       int lab0 = mlt->GetLabel(itr,0);
1091       int lab1 = mlt->GetLabel(itr,1);
1092       int typeMC = 2; // comb.bg.
1093       if (lab0 == lab1) typeMC = fStack->IsPhysicalPrimary(lab0) ? 0:1; // prim or sec
1094       if      (typeMC==0) FillHistosSet(fHistosTrPrim,eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist); // primary
1095       else if (typeMC==1) FillHistosSet(fHistosTrSec, eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist); // secondary
1096       else {
1097         FillHistosSet(fHistosTrComb,eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist); // comb
1098         // for combinatorals fill also the uncorrelated part
1099         if (fMultReco) {
1100           float *trl = fMultReco->GetTracklet(itr);
1101           int clId0 = (int)trl[AliITSMultReconstructor::kClID1];
1102           int clId1 = (int)trl[AliITSMultReconstructor::kClID2];
1103           float *clLabs0 = fMultReco->GetClusterOfLayer(0,clId0) + AliITSMultReconstructor::kClMC0;
1104           float *clLabs1 = fMultReco->GetClusterOfLayer(1,clId1) + AliITSMultReconstructor::kClMC0;
1105           if (!HaveCommonParent(clLabs0,clLabs1))
1106             FillHistosSet(fHistosTrCombU,eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist);
1107         }
1108       } // combinatorials
1109
1110       if (dist<fNStdCut) {
1111         double dphiS  = TMath::Abs(dphi) - fDPhiShift; if (dphi<0) dphiS = -dphiS;
1112         if (dphiS<fDPhiSCut) FillSpecies(typeMC, lab0);
1113       }
1114       if (fCheckReconstructables) CheckReconstructables();
1115     }
1116   }
1117   //
1118   //-------------------------------------------------------------TMP RS - singles ------->>>
1119   int offsH = fCurrCentBin*kNStandardH;
1120   TH2* hSingles = (TH2*)histos->UncheckedAt(offsH+kHEtaZvSPD1);
1121   if (hSingles) {
1122     int nclS = mlt->GetNumberOfSingleClusters();
1123     double *thtS = mlt->GetThetaSingle();
1124     for (int ics=nclS;ics--;) {
1125       double etaS = -TMath::Log(TMath::Tan(thtS[ics]/2));
1126       if (etaS<fEtaMin || etaS>fEtaMax) continue;
1127       hSingles->Fill(etaS,fESDVtx[2]);
1128     }
1129   }
1130   //-------------------------------------------------------------TMP RS - singles -------<<<
1131   //
1132 }
1133
1134 //_________________________________________________________________________
1135 void AliTrackletTaskMulti::FillMCPrimaries()
1136 {
1137   // fill all MC primaries Zv vs Eta
1138   if (!fStack || !fMCEvent) return;
1139
1140   //---------------------------------------- CHECK ------------------------------>>>
1141   TArrayF vtxMC;
1142   AliGenHijingEventHeader* pyHeader = 0;
1143   //
1144   if (fUseMC) {
1145     pyHeader = (AliGenHijingEventHeader*) fMCEvent->GenEventHeader();//header->GenEventHeader();
1146     pyHeader->PrimaryVertex(vtxMC);
1147   }
1148   //---------------------------------------- CHECK ------------------------------<<<
1149   //
1150   int ntr = fStack->GetNtrack();
1151   TH2* hprimEtaZ = (TH2F*)fHistosCustom->UncheckedAt(kHZVEtaPrimMC+fCurrCentBin);
1152   int nprim = 0;
1153   for (int itr=ntr;itr--;) {
1154     if (!fStack->IsPhysicalPrimary(itr)) continue;
1155     AliMCParticle *part  = (AliMCParticle*)fMCEvent->GetTrack(itr);
1156     if (!part->Charge()) continue;
1157     //
1158     //---------------------------------------- CHECK ------------------------------>>>
1159     /*
1160     Float_t dz = part->Zv() - vtxMC[2];
1161     if (TMath::Abs(dz)>1e-6) continue; // reject
1162     */
1163     //---------------------------------------- CHECK ------------------------------<<<
1164     //
1165     Float_t theta = part->Theta();
1166     if (theta<1e-6 || theta>TMath::Pi()-1e-6) continue;
1167     Float_t eta = part->Eta();
1168     if (eta<fEtaMin || eta>fEtaMax) continue;
1169     hprimEtaZ->Fill(eta, fESDVtx[2]);
1170     nprim++;
1171   }
1172   //
1173   ((TH1*)fHistosCustom->UncheckedAt(kHNPrimMeanMC))->Fill(fCurrCentBin,nprim);
1174   if (fNPart>0) {
1175     ((TH1*)fHistosCustom->UncheckedAt(kHNPrim2PartMC))->Fill(fCurrCentBin,nprim/fNPart);
1176     ((TH2*)fHistosCustom->UncheckedAt(kHNPrim2PartNpMC))->Fill(fNPart,nprim/fNPart);
1177     ((TH2*)fHistosCustom->UncheckedAt(kHNPartMC))->Fill(fNPart,fCurrCentBin);
1178     ((TH1*)fHistosCustom->UncheckedAt(kHNPartMeanMC))->Fill(fCurrCentBin,fNPart);
1179   }
1180   if (fNBColl>0) {
1181     ((TH1*)fHistosCustom->UncheckedAt(kHNPrim2BCollMC))->Fill(fCurrCentBin,nprim/fNBColl);
1182     ((TH2*)fHistosCustom->UncheckedAt(kHNPrim2BCollNpMC))->Fill(fNPart,nprim/fNBColl);
1183     ((TH2*)fHistosCustom->UncheckedAt(kHNBCollMC))->Fill(fNBColl,fCurrCentBin);
1184     ((TH1*)fHistosCustom->UncheckedAt(kHNBCollMeanMC))->Fill(fCurrCentBin,fNBColl);
1185   }
1186   //
1187 }
1188
1189 //_________________________________________________________________________
1190  void AliTrackletTaskMulti::FillHistosSet(TObjArray* histos, double eta,
1191                                           //double /*phi*/,double /*theta*/,
1192                                           double dphi,double dtheta,double dThetaX,
1193                                           double dist)
1194 {
1195   // fill standard set of histos
1196   if (dist>fNStdDev) return;
1197   //
1198   double dphiS  = TMath::Abs(dphi) - fDPhiShift;
1199   if (dphi<0) dphiS = -dphiS;
1200   //
1201   int offs = fCurrCentBin*kNStandardH;
1202   //
1203   if (histos->UncheckedAt(offs+kHDPhiSDThetaX))
1204     ((TH2*)histos->UncheckedAt(offs+kHDPhiSDThetaX))->Fill(dphiS,dThetaX);
1205   //
1206   if (histos->UncheckedAt(offs+kHDPhiDTheta))
1207     ((TH2*)histos->UncheckedAt(offs+kHDPhiDTheta))->Fill(dphi,dtheta);
1208   //
1209   if (histos->UncheckedAt(kHWDist))
1210     ((TH2*)histos->UncheckedAt(offs+kHWDist))->Fill(dist);
1211   //
1212   if (dist<fNStdCut && dphiS<fDPhiSCut && histos->UncheckedAt(offs+kHEtaZvCut))
1213     ((TH2*)histos->UncheckedAt(offs+kHEtaZvCut))->Fill(eta,fESDVtx[2]);
1214   //
1215 }
1216
1217 //__________________________________________________________________
1218 void AliTrackletTaskMulti::FillSpecies(Int_t primsec, Int_t id)
1219 {
1220   // fill PDGcode
1221   TH1 *hPart=0,*hParent=0;
1222   if (primsec==0) {
1223     hPart   = (TH1*)fHistosCustom->UncheckedAt(kHPrimPDG);
1224     hParent = (TH1*)fHistosCustom->UncheckedAt(kHPrimParPDG);
1225   }
1226   else if (primsec==1) {
1227     hPart   = (TH1*)fHistosCustom->UncheckedAt(kHSecPDG);
1228     hParent = (TH1*)fHistosCustom->UncheckedAt(kHSecParPDG);
1229   }
1230   else return;
1231   int ntr = fStack->GetNtrack();
1232   TParticle* part = fStack->Particle(id);
1233   int pdgCode = TMath::Abs(part->GetPdgCode());
1234   int pdgBin = GetPdgBin(pdgCode);
1235   int parID = part->GetFirstMother();
1236   int pdgCodePar = -1;
1237   int pdgBinPar = -1;
1238   while (parID>=0 && parID<ntr) {
1239     part = fStack->Particle(parID);
1240     pdgCodePar = TMath::Abs(part->GetPdgCode());
1241     parID = part->GetFirstMother();
1242   }
1243   if (pdgCodePar>0) pdgBinPar = GetPdgBin(pdgCodePar);
1244   //
1245   hPart->Fill(pdgBin,fCurrCentBin);
1246   hParent->Fill(pdgBinPar,fCurrCentBin);
1247   //
1248 }
1249
1250 //_________________________________________________________________________
1251 Int_t AliTrackletTaskMulti::GetCentralityBin(Float_t perc) const
1252 {
1253   // calculate centrality bin
1254   for (int i=0;i<fNCentBins;i++) if (perc>=fgkCentPerc[i] && perc<=fgkCentPerc[i+1]) return i;
1255   return -1;
1256 }
1257
1258 //_________________________________________________________________________
1259 Int_t AliTrackletTaskMulti::GetPdgBin(Int_t pdgCode)
1260 {
1261   // return my pdg bin
1262   int ncodes = sizeof(fgkPDGCodes)/sizeof(int);
1263   int pdgBin=0;
1264   for (pdgBin=0;pdgBin<ncodes;pdgBin++) if (pdgCode==fgkPDGCodes[pdgBin]) break;
1265   if (pdgBin>=ncodes) {
1266     if (float(pdgCode)>1e9) pdgBin = ncodes; // nuclei
1267     else pdgBin = ncodes+1; // unknown
1268   }
1269   return pdgBin;
1270 }
1271
1272 //_________________________________________________________________________
1273 Bool_t AliTrackletTaskMulti::HaveCommonParent(const float* clLabs0,const float* clLabs1)
1274 {
1275   // do 2 clusters have common parrent
1276   const int kMaxPar = 50;
1277   static int pars[2][50];
1278   int npars[2]={0,0};
1279   const float *labs[2] = {clLabs0,clLabs1};
1280   int ntr = fStack->GetNtrack();
1281   for (int il=0;il<2;il++) {
1282     for (int ilb=0;ilb<3;ilb++) {
1283       int lbl = (int)labs[il][ilb];
1284       if (lbl<0 || lbl>=ntr) continue;
1285       //
1286       while (npars[il]<kMaxPar-1) {
1287         pars[il][ npars[il]++ ] = lbl;
1288         TParticle* part = fStack->Particle(lbl);
1289         if (!part) break;
1290         lbl = part->GetFirstMother();
1291         if (lbl<1 || lbl>=ntr) break;
1292       }
1293     }
1294   }
1295   // compare array of parents
1296   for (int i0=npars[0];i0--;) for (int i1=npars[1];i1--;) if (pars[0][i0]==pars[1][i1]) return kTRUE;
1297   return kFALSE;
1298 }
1299
1300
1301 //_________________________________________________________________________
1302 void AliTrackletTaskMulti::CheckReconstructables()
1303 {
1304   // fill reconstructable tracklets hitsos
1305   static TArrayI trInd;
1306   static TBits   isPrimArr;
1307   //
1308   if (!fMultReco || !fMultReco->IsRecoDone()) {AliInfo("To check reconstructables the reco had to be requested"); return;}
1309   if (!fStack) {AliInfo("MC Stack is not availalble"); return;}
1310   const double kPtMin = 0.05;
1311   //
1312   TClonesArray *clArr[2];
1313   for (int ilr=0;ilr<2;ilr++) {
1314     clArr[ilr] = fMultReco->GetClustersOfLayer(ilr);
1315     if (!clArr[ilr]) {AliInfo("Clusters are not available"); return;}
1316   }
1317   //
1318   int ntr = fStack->GetNtrack();
1319   if (!ntr) return;
1320   trInd.Reset();
1321   if (trInd.GetSize()<ntr) trInd.Set(ntr);
1322   isPrimArr.ResetAllBits();
1323   // count track wich may be reconstructable
1324   //
1325   int ntrStore=0,ntrStorePrim=0;
1326   Int_t *trIndArr = trInd.GetArray();
1327   for (Int_t it=0; it<ntr; it++) {
1328     TParticle* part = fStack->Particle(it);
1329     if (TMath::Abs(part->Eta())>2.2)       continue;
1330     if (TMath::Abs(part->Pt())<kPtMin)      continue;
1331     if (fStack->IsPhysicalPrimary(it)) {
1332       isPrimArr.SetBitNumber(it);
1333       ntrStorePrim++;
1334     }
1335     else { // check if secondary is worth cheking
1336       TParticlePDG* pdgPart = part->GetPDG();
1337       if (TMath::Abs(pdgPart->Charge())!=3)  continue;
1338       if (part->R()>5.)                      continue;
1339     }
1340     trIndArr[it] = ++ntrStore;
1341   }
1342   //
1343   AliInfo(Form("Selected %d MC particles (%d prim) out of %d in the stack\n",ntrStore,ntrStorePrim,ntr));
1344   //
1345   const int kMaxCl=3;
1346   AliITSRecPoint **clIndL[2];
1347   clIndL[0] = new AliITSRecPoint*[kMaxCl*ntrStore]; // max 2 clusters per layer
1348   clIndL[1] = new AliITSRecPoint*[kMaxCl*ntrStore]; // max 2 clusters per layer
1349   memset(clIndL[0],0,kMaxCl*ntrStore*sizeof(AliITSRecPoint*));
1350   memset(clIndL[1],0,kMaxCl*ntrStore*sizeof(AliITSRecPoint*));
1351   //
1352   for (int ilr=0;ilr<2;ilr++) {
1353     TClonesArray *clusters = clArr[ilr];
1354     int ncl = clusters->GetEntriesFast();
1355     for (int icl=ncl;icl--;) {
1356       AliITSRecPoint *cl = (AliITSRecPoint*)clusters->UncheckedAt(icl);
1357       for (int ilb=3;ilb--;) {
1358         int lbl = cl->GetLabel(ilb); if (lbl<0 || lbl>=ntr) continue;
1359         int lblI = trIndArr[lbl];
1360         if (--lblI<0) continue;    // not kept
1361         for (int icc=0;icc<kMaxCl;icc++) if (!clIndL[ilr][lblI+icc*ntrStore]) {clIndL[ilr][lblI+ntrStore*icc] = cl; break;} // first empty one
1362       }
1363     }
1364   }
1365   //
1366   Float_t clusterLay[2][AliITSMultReconstructor::kClNPar];
1367   double trComp[6][kMaxCl*kMaxCl];
1368   int indQual[kMaxCl*kMaxCl];
1369   //
1370   for (int itr=ntr;itr--;) {
1371     int lblI = trIndArr[itr];
1372     if (--lblI<0) continue; // discarded
1373     int ntrCand = 0;        // number of tracklet candidates (including overlaps)
1374     for (int icl0=0;icl0<kMaxCl;icl0++) {
1375       AliITSRecPoint *cl0 = clIndL[0][lblI+icl0*ntrStore];
1376       if (!cl0 || !clIndL[1][lblI]) break;
1377       cl0->GetGlobalXYZ( clusterLay[0] );
1378       fMultReco->ClusterPos2Angles(clusterLay[0], fESDVtx);
1379       for (int icl1=0;icl1<kMaxCl;icl1++) {
1380         AliITSRecPoint *cl1 = clIndL[1][lblI+icl1*ntrStore];
1381         if (!cl1) break;
1382         cl1->GetGlobalXYZ( clusterLay[1] );
1383         fMultReco->ClusterPos2Angles(clusterLay[1], fESDVtx);
1384         trComp[AliITSMultReconstructor::kTrPhi][ntrCand]    = clusterLay[0][AliITSMultReconstructor::kClPh];
1385         trComp[AliITSMultReconstructor::kTrTheta][ntrCand]  = clusterLay[0][AliITSMultReconstructor::kClTh];
1386         trComp[AliITSMultReconstructor::kTrDTheta][ntrCand] = clusterLay[0][AliITSMultReconstructor::kClTh] - clusterLay[1][AliITSMultReconstructor::kClTh];
1387         trComp[AliITSMultReconstructor::kTrDPhi][ntrCand]   = clusterLay[0][AliITSMultReconstructor::kClPh] - clusterLay[1][AliITSMultReconstructor::kClPh];
1388         trComp[AliITSMultReconstructor::kTrLab1][ntrCand]   = icl1*10 + icl0;
1389         double &dphi = trComp[ntrCand][3];
1390         if (dphi>TMath::Pi()) dphi=2.*TMath::Pi()-dphi;     // take into account boundary condition
1391         trComp[5][ntrCand] = fMultReco->CalcDist(trComp[AliITSMultReconstructor::kTrDPhi][ntrCand],
1392                                                  trComp[AliITSMultReconstructor::kTrDTheta][ntrCand],
1393                                                  trComp[AliITSMultReconstructor::kTrTheta][ntrCand]);
1394         ntrCand++;
1395       }
1396     }
1397     if (!ntrCand) continue; // no tracklets
1398     if (ntrCand>1) TMath::Sort(ntrCand,trComp[5],indQual,kFALSE); else indQual[0] = 0; // sort in weighted distance
1399     if (fRemoveOverlaps) ntrCand = 1; // select the best
1400     //
1401     // disable worst tracklet with shared cluster
1402     for (int itc=0;itc<ntrCand;itc++) {
1403       int ind = indQual[itc];
1404       if (trComp[AliITSMultReconstructor::kTrLab1][ind]<0) continue; // already disabled
1405       for (int jtc=itc+1;jtc<ntrCand;jtc++) {
1406         int jnd = indQual[jtc];
1407         if (trComp[AliITSMultReconstructor::kTrLab1][jnd]<0) continue; // already disabled
1408         if ( int(trComp[AliITSMultReconstructor::kTrLab1][ind])/10 == int(trComp[AliITSMultReconstructor::kTrLab1][jnd])/10 ||
1409              int(trComp[AliITSMultReconstructor::kTrLab1][ind])%10 == int(trComp[AliITSMultReconstructor::kTrLab1][jnd])%10) trComp[AliITSMultReconstructor::kTrLab1][jnd] = -1;
1410       }
1411     }
1412     //
1413     // store, but forbid cluster reusing
1414     TObjArray* histos = isPrimArr.TestBitNumber(itr) ? fHistosTrRcblPrim : fHistosTrRcblSec;
1415     for (int itc=0;itc<ntrCand;itc++) {
1416       int ind = indQual[itc];
1417       if (trComp[4][ind]<0) continue; // discarded
1418       double eta    = -TMath::Log(TMath::Tan(trComp[AliITSMultReconstructor::kTrTheta][ind]/2));
1419       if (eta<fEtaMin || eta>fEtaMax) continue;
1420       double dThetaX = trComp[AliITSMultReconstructor::kTrTheta][ind];
1421       if (fScaleDTBySin2T) {
1422         double sint   =  TMath::Sin(trComp[AliITSMultReconstructor::kTrTheta][ind]);
1423         dThetaX /= (sint*sint);
1424       }
1425       FillHistosSet(histos,eta,
1426                     //trComp[AliITSMultReconstructor::kTrPhi][ind],trComp[AliITSMultReconstructor::kTrTheta][ind],
1427                     trComp[AliITSMultReconstructor::kTrDPhi][ind],trComp[AliITSMultReconstructor::kTrDTheta][ind],
1428                     dThetaX,trComp[5][ind]);
1429     }
1430   }
1431   //
1432   delete[] clIndL[0];
1433   delete[] clIndL[1];
1434 }
1435
1436 //_________________________________________________________________________
1437 void AliTrackletTaskMulti::FillClusterInfo()
1438 {
1439   // fill info on clusters associated to good tracklets
1440   if (!fMultReco) return;
1441   int ntr = fMultReco->GetNTracklets();
1442   int clID[2];
1443   TH2F *hclU[2] = {(TH2F*)fHistosCustom->UncheckedAt(kHClUsedInfoL0),(TH2F*)fHistosCustom->UncheckedAt(kHClUsedInfoL1)};
1444   TH2F *hclA[2] = {(TH2F*)fHistosCustom->UncheckedAt(kHClAllInfoL0),(TH2F*)fHistosCustom->UncheckedAt(kHClAllInfoL1)};
1445   for (int itr=ntr;itr--;) {
1446     Float_t *trc = fMultReco->GetTracklet(itr);
1447     if (TMath::Abs(TMath::Abs(trc[AliITSMultReconstructor::kTrDPhi])-fDPhiShift)>fDPhiSCut) continue;
1448     if (fMultReco->CalcDist(trc[AliITSMultReconstructor::kTrDPhi],
1449                             trc[AliITSMultReconstructor::kTrDTheta],
1450                             trc[AliITSMultReconstructor::kTrTheta]) > fNStdCut) continue;
1451     clID[0] = (int)trc[AliITSMultReconstructor::kClID1];
1452     clID[1] = (int)trc[AliITSMultReconstructor::kClID2];
1453     for (int il=0;il<2;il++) {
1454       Float_t *clinf = fMultReco->GetClusterOfLayer(il,clID[il]);
1455       hclU[il]->Fill( clinf[AliITSMultReconstructor::kClZ], clinf[AliITSMultReconstructor::kClPh]);
1456     }
1457   }
1458   //
1459   for (int il=0;il<2;il++) for (int ic=fMultReco->GetNClustersLayer(il);ic--;) {
1460       Float_t *clinf = fMultReco->GetClusterOfLayer(il,ic);
1461       hclA[il]->Fill( clinf[AliITSMultReconstructor::kClZ], clinf[AliITSMultReconstructor::kClPh]);
1462     }
1463   //
1464 }
1465
1466 //_________________________________________________________________________
1467 void AliTrackletTaskMulti::FillClusterInfoFromMult(const AliMultiplicity* mlt, double zVertex)
1468 {
1469   // fill info on clusters taking them from Multiplicity object
1470   const double kRSPD2 = 3.9;
1471   TH2F *hclU = (TH2F*)fHistosCustom->UncheckedAt(kHClUsedInfoL0);
1472   TH2F *hclA = (TH2F*)fHistosCustom->UncheckedAt(kHClAllInfoL0);
1473   int ntr = mlt->GetNumberOfTracklets();
1474   for (int itr=ntr;itr--;) {
1475     Bool_t goodTracklet = kTRUE;
1476     if (TMath::Abs( mlt->GetDeltaPhi(itr)-fDPhiShift)>fDPhiSCut) goodTracklet = kFALSE;
1477     if (mlt->CalcDist(itr) > fNStdCut) goodTracklet = kFALSE;
1478     double phi   = mlt->GetPhi(itr);
1479     double z     = kRSPD2/TMath::Tan(mlt->GetTheta(itr)) + zVertex;
1480     if (goodTracklet) hclU->Fill(z,phi);
1481     hclA->Fill(z,phi);
1482   }
1483   //
1484   int ncl = mlt->GetNumberOfSingleClusters();
1485   for (int icl=ncl;icl--;) {
1486     double phi   = mlt->GetPhiSingle(icl);
1487     double z     = kRSPD2/TMath::Tan(mlt->GetThetaSingle(icl)) + zVertex;
1488     hclA->Fill(z,phi);
1489   }
1490   //
1491 }
1492