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