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