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
  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::fgCentSelName[] = {"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 (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) {  //RRR
 277     printf("Deleteing output\n");
 278     delete fOutput;
 279     fOutput = 0;
 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 AliTrackletTaskMulti::UserCreateOutputObjects()
 300 {
 301   //
 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",fgCentSelName[fUseCentralityVar]));
 371   PostData(1, fOutput);
 372   //
 373 }
 374
 375 //________________________________________________________________________
 376 void AliTrackletTaskMulti::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(fgCentSelName[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   }
 539   //
 540   // Injection: it must come right after the normal reco since needs its results
 541   if (GetDoInjection()) {
 542     if (!fMultReco) InitMultReco(); // in principle, not needed, the reco is created above
 543     fMultReco->SetRecType(AliITSMultRecBg::kBgInj);
 544     fMultReco->Run(fRPTree, vtxf);
 545     printf("Multiplicity from Injection:\n");
 546     AliMultiplicity* mlt = fMultReco->GetMultiplicity();
 547     if (mlt) mlt->Print();
 548     hstat->Fill(kBinEntries + kEvProcInj + kEntriesPerBin*fCurrCentBin);
 549     FillHistos(kBgInj,mlt);
 550   }
 551   //
 552   // Rotation
 553   if (GetDoRotation()) {
 554     InitMultReco();
 555     fMultReco->SetRecType(AliITSMultRecBg::kBgRot);
 556     fMultReco->SetPhiRotationAngle(fPhiRot);
 557     fMultReco->Run(fRPTree, vtxf);
 558     printf("Multiplicity from Rotation:\n");
 559     AliMultiplicity* mlt = fMultReco->GetMultiplicity();
 560     if (mlt) mlt->Print();
 561     hstat->Fill(kBinEntries + kEvProcRot + kEntriesPerBin*fCurrCentBin);
 562     FillHistos(kBgRot,mlt);
 563   }
 564   //
 565   if (GetDoMixing()) {
 566     /*
 567     AliMixEventInputHandler* handToMix = (AliMixEventInputHandler*)handRP->MixingHandler();
 568     if (!handToMix) { printf("No Mixing handler\n"); return; }
 569     handToMix->GetEntry();
 570     if(handToMix->MixedEventNumber()<1) {printf("Mixing: No enough events in pool\n"); return;}
 571     AliESDInputHandlerRP* handRPMix = (AliESDInputHandlerRP*) handToMix->InputEventHandler(0);
 572
 573     if (!handRPMix) { printf("No Mixing RP handler\n"); return; }
 574     fRPTreeMix = handRPMix->GetTreeR("ITS");
 575     if (!fRPTreeMix) { AliError(" Invalid ITS cluster tree of the 2nd event!\n"); return; }
 576     //
 577     AliESDEvent *esdMix = handRPMix->GetEvent();
 578     const AliESDVertex* vtxESDmix = esdMix->GetVertex();
 579     ((TH2*)fHistosCustom->UncheckedAt(kHZVtxMixDiff))->Fill(vtxESDmix->GetZ()-esdvtx[2],fCurrCentBin);
 580     ((TH2*)fHistosCustom->UncheckedAt(kHNTrMixDiff) )->
 581       Fill(esdMix->GetMultiplicity()->GetNumberOfTracklets() - multESD->GetNumberOfTracklets(),fCurrCentBin);
 582     //
 583     InitMultReco();
 584     fMultReco->SetRecType(AliITSMultRecBg::kBgMix);
 585     fMultReco->Run(fRPTree, vtxf,fRPTreeMix);
 586     printf("Multiplicity from Mixing:\n");
 587     AliMultiplicity* mlt = fMultReco->GetMultiplicity();
 588     if (mlt) mlt->Print();
 589     hstat->Fill(kBinEntries + kEvProcMix + kEntriesPerBin*fCurrCentBin);
 590     FillHistos(kBgMix,mlt);
 591     */
 592     //
 593   }
 594   // =============================================================================<<<
 595   //
 596   if (fMultReco) delete fMultReco;
 597   fMultReco = 0;
 598   //
 599 }
 600
 601
 602 //________________________________________________________________________
 603 void AliTrackletTaskMulti::Terminate(Option_t *)
 604 {
 605   Printf("Terminating...");
 606   TH1* hstat;
 607   TList *lst = dynamic_cast<TList*>(GetOutputData(1));
 608   printf("Term: %p %p %p\n",fOutput,lst,fHistosCustom);
 609   if (lst && (hstat=(TH1*)lst->FindObject("hStat"))) {
 610     Info("Terminate","Registering used settings");
 611     // fill used settings
 612     hstat->Fill(kOneUnit,1.);
 613     hstat->Fill(kCentVar,fUseCentralityVar);
 614     hstat->Fill(kDPhi,fDPhiWindow);
 615     hstat->Fill(kDTht,fDThetaWindow);
 616     hstat->Fill(kNStd,fNStdDev);
 617     hstat->Fill(kPhiShift,fDPhiShift);
 618     hstat->Fill(kThtS2,fScaleDTBySin2T);
 619     hstat->Fill(kThtCW,fCutOnDThetaX);
 620     hstat->Fill(kPhiOvl,fPhiOverlapCut);
 621     hstat->Fill(kZEtaOvl,fZetaOverlap);
 622     hstat->Fill(kNoOvl,fRemoveOverlaps);
 623     //
 624     hstat->Fill(kPhiRot,fPhiRot);
 625     hstat->Fill(kInjScl,fInjScale);
 626     hstat->Fill(kEtaMin,fEtaMin);
 627     hstat->Fill(kEtaMax,fEtaMax);
 628     hstat->Fill(kZVMin,fZVertexMin);
 629     hstat->Fill(kZVMax,fZVertexMax);
 630     //
 631     hstat->Fill(kDPiSCut,fDPhiSCut);
 632     hstat->Fill(kNStdCut,fNStdCut);
 633     hstat->Fill(kMCV0Scale, fMCV0Scale);
 634     //
 635   }
 636   //
 637   //  AliAnalysisTaskSE::Terminate();
 638 }
 639
 640
 641 //_________________________________________________________________________
 642 void AliTrackletTaskMulti::InitMultReco()
 643 {
 644   // create mult reconstructor
 645   if (fMultReco) delete fMultReco;
 646   fMultReco = new AliITSMultRecBg();
 647   fMultReco->SetCreateClustersCopy(kTRUE);
 648   fMultReco->SetScaleDThetaBySin2T(fScaleDTBySin2T);
 649   fMultReco->SetNStdDev(fNStdDev);
 650   fMultReco->SetPhiWindow( fDPhiWindow );
 651   fMultReco->SetThetaWindow( fDThetaWindow );
 652   fMultReco->SetPhiShift( fDPhiShift );
 653   fMultReco->SetRemoveClustersFromOverlaps(fRemoveOverlaps);
 654   fMultReco->SetPhiOverlapCut(fPhiOverlapCut);
 655   fMultReco->SetZetaOverlapCut(fZetaOverlap);
 656   fMultReco->SetHistOn(kFALSE);
 657   fMultReco->SetRecType( AliITSMultRecBg::kData );
 658 }
 659
 660 //_________________________________________________________________________
 661 TObjArray* AliTrackletTaskMulti::BookCustomHistos()
 662 {
 663   // book custom histos, not related to specific tracklet type
 664   TObjArray* histos = new TObjArray();
 665   TH1F* hstat;
 666   //
 667   // ------------ job parameters, statistics ------------------------------>>>
 668   int nbs = kBinEntries + fNCentBins*kEntriesPerBin;
 669   hstat = new TH1F("hStat","Run statistics",nbs,0.5,nbs+0.5);
 670   //
 671   hstat->GetXaxis()->SetBinLabel(kEvTot, "Ev.Tot0");
 672   hstat->GetXaxis()->SetBinLabel(kEvTot, "Ev.Tot");
 673   hstat->GetXaxis()->SetBinLabel(kOneUnit,"ScaleMerge");
 674   hstat->GetXaxis()->SetBinLabel(kNWorkers,"Workers");
 675   //
 676   hstat->GetXaxis()->SetBinLabel(kCentVar,  fgCentSelName[fUseCentralityVar]);
 677   hstat->GetXaxis()->SetBinLabel(kDPhi,  "#Delta#varphi");
 678   hstat->GetXaxis()->SetBinLabel(kDTht,  "#Delta#theta");
 679   hstat->GetXaxis()->SetBinLabel(kNStd,  "N.std");
 680   hstat->GetXaxis()->SetBinLabel(kPhiShift,"#delta#varphi");
 681   hstat->GetXaxis()->SetBinLabel(kThtS2,"scale #Delta#theta");
 682   hstat->GetXaxis()->SetBinLabel(kPhiOvl,"#varpho_{Ovl}");
 683   hstat->GetXaxis()->SetBinLabel(kZEtaOvl,"#z_{Ovl}");
 684   hstat->GetXaxis()->SetBinLabel(kNoOvl, "rem.ovl");
 685   //
 686   hstat->GetXaxis()->SetBinLabel(kPhiRot,"#varphi_{rot}");
 687   hstat->GetXaxis()->SetBinLabel(kInjScl,"inj");
 688   hstat->GetXaxis()->SetBinLabel(kEtaMin,"#eta_{min}");
 689   hstat->GetXaxis()->SetBinLabel(kEtaMax,"#eta_{max}");
 690   hstat->GetXaxis()->SetBinLabel(kZVMin,"ZV_{min} cut");
 691   hstat->GetXaxis()->SetBinLabel(kZVMax,"ZV_{max} cut");
 692   //
 693   hstat->GetXaxis()->SetBinLabel(kDPiSCut,"#Delta#varphi-#delta_{#phi} cut");
 694   hstat->GetXaxis()->SetBinLabel(kNStdCut,"#Delta cut");
 695   //
 696   hstat->GetXaxis()->SetBinLabel(kMCV0Scale,"MC V0 scale");
 697   //
 698   for (int i=0;i<fNCentBins;i++) {
 699     TString bnt = "b"; bnt+= i;
 700     int offs = kBinEntries + i*kEntriesPerBin;
 701     hstat->GetXaxis()->SetBinLabel(offs + kEvProcData, bnt+" Ev.ProcData");
 702     hstat->GetXaxis()->SetBinLabel(offs + kEvProcInj,  bnt+" Ev.ProcInj");
 703     hstat->GetXaxis()->SetBinLabel(offs + kEvProcRot,  bnt+" Ev.ProcRot");
 704     hstat->GetXaxis()->SetBinLabel(offs + kEvProcMix,  bnt+" Ev.ProcMix");
 705     //
 706   }
 707   //
 708   hstat->Fill(kNWorkers);
 709   //
 710   AddHisto(histos,hstat,kHStat);
 711   //
 712   // ------------------------ events per centrality bin ----------------------
 713   TH1D* hCentAx = new TH1D("EvCentr","Events per centrality",fNCentBins,fgkCentPerc);
 714   hCentAx->GetXaxis()->SetTitle("Centrality parameter");
 715   AddHisto(histos,hCentAx,kHStatCent);
 716   //
 717   TH1D* hCentBin = new TH1D("EvCentrBin","Events per centrality bin",fNCentBins,-0.5,fNCentBins-0.5);
 718   hCentBin->GetXaxis()->SetTitle("Centrality Bin");
 719   AddHisto(histos,hCentBin,kHStatCentBin);
 720   //
 721   // ------------ job parameters, statistics ------------------------------<<<
 722   //
 723   double etaMn=-3,etaMx=3;
 724   double zMn=-30, zMx=30;
 725   int nEtaBins = int((etaMx-etaMn)/0.1);
 726   if (nEtaBins<1) nEtaBins = 1;
 727   //
 728   int nZVBins = int(zMx-zMn);
 729   if (nZVBins<1) nZVBins = 1;
 730   //
 731   // Z vertex distribution for events before selection
 732   TH1F* hzvns = new  TH1F("zvNoSel","Z vertex before selection",nZVBins,zMn,zMx);
 733   hzvns->GetXaxis()->SetTitle("Zvertex");
 734   AddHisto(histos,hzvns,kHZVtxNoSel);
 735   //
 736   // V0 for events before selection
 737   int nbmltV0 = 250;
 738   double maxmltV0 = 25000;
 739   //
 740   TH1F* hnV0ns = new  TH1F("V0NoSel","V0 signal Before Cent. Selection",nbmltV0,0,maxmltV0);
 741   hnV0ns->GetXaxis()->SetTitle("V0 signal");
 742   AddHisto(histos,hnV0ns,kHV0NoSel);
 743   //
 744   // N SPD2 clusters
 745   int nbmltSPD2 = 175;
 746   double maxmltSPD2 = 7000;
 747   TH1F* hncl2ns = new  TH1F("NClustersSPD2NoSel","N Clusters on SPD2 Before Cent Selection",nbmltSPD2,0,maxmltSPD2);
 748   hncl2ns->GetXaxis()->SetTitle("N Clus SPD2");
 749   AddHisto(histos,hncl2ns,kHNClSPD2NoSel);
 750   //
 751   int nbzdc=50;
 752   double maxZDC=6000., maxZEM=2500.;
 753   TH2F* hzdczemns = new  TH2F("ZDCZEMNoSel","ZDC vs ZEM Before Cent Selection",
 754                               nbzdc,0,maxZEM,nbzdc,0,maxZDC);
 755   hzdczemns->GetXaxis()->SetTitle("ZEM");
 756   hzdczemns->GetXaxis()->SetTitle("ZDC");
 757   AddHisto(histos,hzdczemns,kHZDCZEMNoSel);
 758   //
 759   TH2F* hzv = new  TH2F("zv","Z vertex after Selection per Cent.Bin",nZVBins,zMn,zMx, fNCentBins, -0.5,fNCentBins-0.5);
 760   hzv->GetXaxis()->SetTitle("Zvertex");
 761   hzv->GetYaxis()->SetTitle("Cent.Bin ID");
 762   AddHisto(histos,hzv,kHZVtx);
 763   //
 764   // V0
 765   TH2F* hnV0 = new  TH2F("V0","V0 signal per Cent.Bin ",nbmltV0,0,maxmltV0, fNCentBins, -0.5,fNCentBins-0.5);
 766   hnV0->GetXaxis()->SetTitle("V0 signal");
 767   hnV0->GetYaxis()->SetTitle("Cent.Bin ID");
 768   AddHisto(histos,hnV0,kHV0);
 769   //
 770   // N SPD2 clusters
 771   TH2F* hncl2 = new  TH2F("NClustersSPD2","N Clusters on SPD2 per Cent.Bin",nbmltSPD2,0,maxmltSPD2, fNCentBins, -0.5,fNCentBins-0.5);
 772   hncl2->GetXaxis()->SetTitle("N Clus SPD2");
 773   hncl2->GetYaxis()->SetTitle("Cent.Bin ID");
 774   AddHisto(histos,hncl2,kHNClSPD2);
 775   //
 776   // ZDCZEM
 777   TH3F* hzdczem = new TH3F("ZDCZEM","ZDC vs ZEM per Cent.Bin",nbzdc,0,maxZEM,nbzdc,0,maxZDC, fNCentBins, -0.5,fNCentBins-0.5);
 778   hzdczem->GetXaxis()->SetTitle("ZEM");
 779   hzdczem->GetYaxis()->SetTitle("ZDC");
 780   AddHisto(histos,hzdczem,kHZDCZEM);
 781   //
 782   //----------------------------------------------------------------------
 783   int nEtaBinsS = int((fEtaMax-fEtaMin)/0.1);
 784   if (nEtaBinsS<1) nEtaBins = 1;
 785   //
 786   int nZVBinsS = int(fZVertexMax-fZVertexMin);
 787   if (nZVBinsS<1) nZVBinsS = 1;
 788
 789   if (fUseMC) {
 790     // Z vertex vs Eta distribution for primaries
 791     char buffn[100],bufft[500];
 792     for (int ib=0;ib<fNCentBins;ib++) {
 793       sprintf(buffn,"b%d_zvEtaPrimMC",ib);
 794       sprintf(bufft,"bin%d Zvertex vs #eta PrimMC",ib);
 795       TH2F* hzvetap = new  TH2F(buffn,bufft, nEtaBinsS,fEtaMin,fEtaMax,nZVBinsS,fZVertexMin,fZVertexMax);
 796       hzvetap->GetXaxis()->SetTitle("#eta");
 797       hzvetap->GetYaxis()->SetTitle("Zvertex");
 798       AddHisto(histos,hzvetap,kHZVEtaPrimMC+ib);
 799     }
 800     //
 801     // <n> primaries according to MC generator
 802     TH1F* hnprimM = new  TH1F("nPrimMean","<N> primaries",fNCentBins, -0.5,fNCentBins-0.5);
 803     hnprimM->GetXaxis()->SetTitle("Cent.Bin ID");
 804     AddHisto(histos,hnprimM,kHNPrimMeanMC);
 805     //
 806     // <n> primaries per part.pair according to MC generator
 807     TH1F* hnprim2partM = new  TH1F("nPrim2Part","<N> primaries per part.pair",fNCentBins, -0.5,fNCentBins-0.5);
 808     hnprim2partM->GetXaxis()->SetTitle("Cent.Bin ID");
 809     AddHisto(histos,hnprim2partM,kHNPrim2PartMC);
 810     //
 811     // <n> primaries per part.pair vs npart.pair according to MC generator
 812     TH2F* hnprim2partNp = new  TH2F("nPrim2Part_vs_NPart","<N> primaries per part.pair vs N part.pairs",105,0,210,200,0,40);
 813     hnprim2partNp->GetXaxis()->SetTitle("N.part.pairs");
 814     hnprim2partNp->GetYaxis()->SetTitle("N.prim/N.part.pairs");
 815     AddHisto(histos,hnprim2partNp,kHNPrim2PartNpMC);
 816     //
 817     // <n> primaries per b.coll vs npart.pair according to MC generator
 818     TH2F* hnprim2BCollNp = new  TH2F("nPrim2BColl_vs_NPart","<N> primaries per bin.coll vs N part.pairs",105,0,210,200,0,40);
 819     hnprim2BCollNp->GetXaxis()->SetTitle("N.part.pairs");
 820     hnprim2BCollNp->GetYaxis()->SetTitle("N.prim/N.bin.coll.");
 821     AddHisto(histos,hnprim2BCollNp,kHNPrim2BCollNpMC);
 822     //
 823     // <n> primaries per bin.coll. according to MC generator
 824     TH1F* hnprim2BCollM = new  TH1F("nPrim2BColl","<N> primaries per bin.coll",fNCentBins, -0.5,fNCentBins-0.5);
 825     hnprim2BCollM->GetXaxis()->SetTitle("Cent.Bin ID");
 826     AddHisto(histos,hnprim2BCollM,kHNPrim2BCollMC);
 827     //
 828     // n participants according to MC generator
 829     TH2F* hnpart = new  TH2F("nPart","N participant pairs",210,0,210,fNCentBins, -0.5,fNCentBins-0.5);
 830     hnpart->GetXaxis()->SetTitle("N part. pairs");
 831     hnpart->GetYaxis()->SetTitle("Cent.Bin ID");
 832     AddHisto(histos,hnpart,kHNPartMC);
 833     //
 834     // <n> participants according to MC generator
 835     TH1F* hnpartM = new  TH1F("nPartMean","<N> participant pairs",fNCentBins, -0.5,fNCentBins-0.5);
 836     hnpartM->GetXaxis()->SetTitle("Cent.Bin ID");
 837     AddHisto(histos,hnpartM,kHNPartMeanMC);
 838     //
 839     // n bin coll. according to MC generator
 840     TH2F* hnbcoll = new  TH2F("nBColl","N bin. coll",2000,0,2000,fNCentBins, -0.5,fNCentBins-0.5);
 841     hnbcoll->GetXaxis()->SetTitle("N bin. coll");
 842     hnbcoll->GetYaxis()->SetTitle("Cent.Bin ID");
 843     AddHisto(histos,hnbcoll,kHNBCollMC);
 844     //
 845     // <n> bin col according to MC generator
 846     TH1F* hnbcollM = new  TH1F("nBCollMean","<N> bin.colls",fNCentBins, -0.5,fNCentBins-0.5);
 847     hnbcollM->GetXaxis()->SetTitle("Cent.Bin ID");
 848     AddHisto(histos,hnbcollM,kHNBCollMeanMC);
 849     //
 850   }
 851   //
 852   if (GetDoMixing()) {
 853     //
 854     // Difference in Z vertex for mixed events
 855     TH2F* hzdiff = new TH2F("MixSPDVertexDiff","SPD #Delta Z Vertex distribution per mult bin ",100,-5,5, fNCentBins, -0.5,fNCentBins-0.5);
 856     hzdiff->GetXaxis()->SetTitle("#Delta Z Vertex [cm]");
 857     hzdiff->GetYaxis()->SetTitle(Form("Entries / %1.2f [cm] per mult bin",10./100.));
 858     AddHisto(histos,hzdiff,kHZVtxMixDiff);
 859     //
 860     // Difference in N tracklets for mixed events
 861     TH2F* hntdiff = new TH2F("MixNTrackletsDiff"," SPD tracklets Diff ",200,-1000,1000, fNCentBins, -0.5,fNCentBins-0.5);
 862     hntdiff->GetXaxis()->SetTitle("# tracklet diff");
 863     AddHisto(histos,hntdiff,kHNTrMixDiff);
 864   }
 865   //
 866   // --------------------------------------------------
 867   if (fUseMC) {
 868     int npdg = sizeof(fgkPDGNames)/sizeof(char*);
 869     TH2F* hpdgP = new TH2F("pdgPrim","primary PDG",npdg,0,npdg,fNCentBins, -0.5,fNCentBins-0.5);
 870     AddHisto(histos,hpdgP,kHPrimPDG);
 871     TH2F* hpdgS = new TH2F("pdgSec","secondary PDG",npdg,0,npdg,fNCentBins, -0.5,fNCentBins-0.5);
 872     AddHisto(histos,hpdgS,kHSecPDG);
 873     TH2F* hpdgPP = new TH2F("pdgPrimPar","primary parent PDG ",npdg,0,npdg,fNCentBins, -0.5,fNCentBins-0.5);
 874     AddHisto(histos,hpdgPP,kHPrimParPDG);
 875     TH2F* hpdgSP = new TH2F("pdgSecPar","secondary parent PDG",npdg,0,npdg,fNCentBins, -0.5,fNCentBins-0.5);
 876     AddHisto(histos,hpdgSP,kHSecParPDG);
 877     for (int i=0;i<npdg;i++) {
 878       hpdgP->GetXaxis()->SetBinLabel(i+1,fgkPDGNames[i]);
 879       hpdgS->GetXaxis()->SetBinLabel(i+1,fgkPDGNames[i]);
 880       hpdgPP->GetXaxis()->SetBinLabel(i+1,fgkPDGNames[i]);
 881       hpdgSP->GetXaxis()->SetBinLabel(i+1,fgkPDGNames[i]);
 882     }
 883   }
 884   //
 885   // -------------------------------------------------
 886   TH2F* hclinf=0;
 887   hclinf = new TH2F("cl0InfoUsed","#phi vs Z of used clusters, Lr0",64,-16,16, 80,0,2*TMath::Pi());
 888   AddHisto(histos,hclinf,kHClUsedInfoL0);
 889   hclinf = new TH2F("cl1InfoUsed","#phi vs Z of used clusters, Lr1",64,-16,16, 2*80,0,2*TMath::Pi());
 890   AddHisto(histos,hclinf,kHClUsedInfoL1);
 891   hclinf = new TH2F("cl0InfoAll","#phi vs Z of all clusters, Lr0",64,-16,16, 80,0,2*TMath::Pi());
 892   AddHisto(histos,hclinf,kHClAllInfoL0);
 893   hclinf = new TH2F("cl1InfoAll","#phi vs Z of all clusters, Lr1",64,-16,16, 2*80,0,2*TMath::Pi());
 894   AddHisto(histos,hclinf,kHClAllInfoL1);
 895   //
 896   // -------------------------------------------------
 897   histos->SetOwner(kFALSE);
 898   //
 899   return histos;
 900 }
 901
 902 //_________________________________________________________________________
 903 TObjArray* AliTrackletTaskMulti::BookHistosSet(const char* pref, UInt_t selHistos)
 904 {
 905   // book standard set of histos attaching the pref in front of the name/title
 906   //
 907   const int kNDPhiBins = 100;
 908   const int kNDThtBins = 100;
 909   int nDistBins = int(fNStdDev)*5;
 910   //
 911   int nEtaBins = int((fEtaMax-fEtaMin)/0.1);
 912   if (nEtaBins<1) nEtaBins = 1;
 913   //
 914   int nZVBins = int(fZVertexMax-fZVertexMin);
 915   if (nZVBins<1) nZVBins = 1;
 916   float dphir = fDPhiWindow*TMath::Sqrt(fNStdDev);
 917   float dthtr = fDThetaWindow*TMath::Sqrt(fNStdDev);
 918   //
 919   TObjArray* histos = new TObjArray();
 920   TH2F* h2;
 921   TH1F* h1;
 922   char buffn[100],bufft[500];
 923   //
 924   for (int ib=0;ib<fNCentBins;ib++) {
 925     //
 926     int offs = ib*kNStandardH;
 927     if (selHistos & (0x1<<kHEtaZvCut) ) {
 928       sprintf(buffn,"b%d_%s_ZvEtaCutT",ib,pref);
 929       sprintf(bufft,"bin%d (%s) Zv vs Eta with tracklet cut",ib,pref);
 930       h2 = new TH2F(buffn,bufft,nEtaBins,fEtaMin,fEtaMax, nZVBins, fZVertexMin,fZVertexMax);
 931       h2->GetXaxis()->SetTitle("#eta");
 932       h2->GetYaxis()->SetTitle("Zv");
 933       AddHisto(histos,h2,offs+kHEtaZvCut);
 934     }
 935     //
 936     if (selHistos & (0x1<<kHDPhiDTheta) ) {
 937       sprintf(buffn,"b%d_%s_dPhidTheta",ib,pref);
 938       sprintf(bufft,"bin%d (%s) #Delta#theta vs #Delta#varphi",ib,pref);
 939       h2 = new TH2F(buffn,bufft,kNDPhiBins,-dphir,dphir,kNDThtBins,-dthtr,dthtr);
 940       h2->GetXaxis()->SetTitle("#Delta#varphi [rad]");
 941       h2->GetYaxis()->SetTitle("#Delta#theta [rad]");
 942       AddHisto(histos,h2,offs+kHDPhiDTheta);
 943     }
 944     //
 945     if (selHistos & (0x1<<kHDPhiSDThetaX) ) {
 946       sprintf(buffn,"b%d_%s_dPhiSdThetaX",ib,pref);
 947       sprintf(bufft,"bin%d (%s) #Delta#theta%s vs #Delta#varphi-#delta_{#varphi}",ib,pref,fScaleDTBySin2T ? "/sin^{2}(#theta)":"");
 948       h2 = new TH2F(buffn,bufft,kNDPhiBins,-dphir,dphir,kNDThtBins,-dthtr,dthtr);
 949       h2->GetXaxis()->SetTitle("#Delta#varphi-#delta_{#varphi} [rad]");
 950       sprintf(bufft,"#Delta#theta%s",fScaleDTBySin2T ? "/sin^{2}(#theta)":"");
 951       h2->GetYaxis()->SetTitle(bufft);
 952       AddHisto(histos,h2,offs+kHDPhiSDThetaX);
 953     }
 954     //
 955     if (selHistos & (0x1<<kHWDist) ) {
 956       sprintf(buffn,"b%d_%s_WDist",ib,pref);
 957       sprintf(bufft,"bin%d #Delta=[(#Delta#varphi-#delta_{#varphi})/#sigma#varphi]^{2}+"
 958               "[#Delta#theta%s/#sigma#theta]^{2}",ib,fScaleDTBySin2T ? "*sin^{-2}(#theta)":"");
 959       h1 = new TH1F(buffn,bufft,nDistBins,0,fNStdDev);
 960       sprintf(bufft,"#Delta=[(#Delta#varphi-#delta_{#varphi})/#sigma#varphi]^{2}+"
 961               "[#Delta#theta%s/#sigma#theta]^{2}",fScaleDTBySin2T ? "*sin^{-2}(#theta)":"");
 962       h1->GetXaxis()->SetTitle(bufft);
 963       AddHisto(histos,h1,offs+kHWDist);
 964     }
 965     //
 966     if (selHistos & (0x1<<kHEtaZvSPD1) ) {
 967       sprintf(buffn,"b%d_%s_ZvEtaSPD1",ib,pref);
 968       sprintf(bufft,"bin%d (%s) Zv vs Eta SPD1 clusters",ib,pref);
 969       h2 = new TH2F(buffn,bufft,nEtaBins,fEtaMin,fEtaMax, nZVBins, fZVertexMin,fZVertexMax);
 970       h2->GetXaxis()->SetTitle("#eta");
 971       h2->GetYaxis()->SetTitle("Zv");
 972       AddHisto(histos,h2,offs+kHEtaZvSPD1);
 973     }
 974     //
 975   }
 976   //
 977   histos->SetOwner(kFALSE);
 978   return histos;
 979 }
 980
 981 //_________________________________________________________________________
 982 void AliTrackletTaskMulti::AddHisto(TObjArray* histos, TObject* h, Int_t at)
 983 {
 984   // add single histo to the set
 985   if (at>=0) histos->AddAtAndExpand(h,at);
 986   else       histos->Add(h);
 987   fOutput->Add(h);
 988 }
 989
 990 //_________________________________________________________________________
 991 void AliTrackletTaskMulti::FillHistos(Int_t type, const AliMultiplicity* mlt)
 992 {
 993   // fill histos of given type
 994   if (!mlt) return;
 995   //
 996   TObjArray* histos = 0;
 997   if      (type == kData)  histos = fHistosTrData;
 998   else if (type == kBgInj) histos = fHistosTrInj;
 999   else if (type == kBgRot) histos = fHistosTrRot;
1000   else if (type == kBgMix) histos = fHistosTrMix;
1001   //
1002   Bool_t fillMC = (type==kData) && fUseMC && fStack;
1003   //
1004   //
1005   //---------------------------------------- CHECK ------------------------------>>>
1006   TArrayF vtxMC;
1007   AliGenHijingEventHeader* pyHeader = 0;
1008   //
1009   if (fUseMC) {
1010     pyHeader = (AliGenHijingEventHeader*) fMCEvent->GenEventHeader();//header->GenEventHeader();
1011     pyHeader->PrimaryVertex(vtxMC);
1012   }
1013   //---------------------------------------- CHECK ------------------------------<<<
1014   //
1015   if (!histos) return;
1016   int ntr = mlt->GetNumberOfTracklets();
1017   for (int itr=ntr;itr--;) {
1018     //
1019     //---------------------------------------- CHECK ------------------------------>>>
1020     /*
1021     if (fUseMC) {
1022       Bool_t reject = kFALSE;
1023       while(1) {
1024         int lab0 = mlt->GetLabel(itr,0);
1025         int lab1 = mlt->GetLabel(itr,1);
1026         if (lab0!=lab1) break;
1027         if (!fStack->IsPhysicalPrimary(lab0)) break;
1028         //
1029         TParticle* part = fStack->Particle(lab0);
1030         Float_t dz = part->Vz() - vtxMC[2];
1031         if (TMath::Abs(dz)<1e-6) break;
1032         reject = kTRUE;
1033         break;
1034       }
1035       if (reject) continue;
1036     }
1037     */
1038     //---------------------------------------- CHECK ------------------------------<<<
1039     //
1040     double theta  = mlt->GetTheta(itr);
1041     double eta    = -TMath::Log(TMath::Tan(theta/2));
1042     if (eta<fEtaMin || eta>fEtaMax) continue;
1043     //
1044     double dtheta = mlt->GetDeltaTheta(itr);
1045     double dThetaX = dtheta;
1046     if (fScaleDTBySin2T) {
1047       double sint   =  TMath::Sin(theta);
1048       dThetaX /= (sint*sint);
1049     }
1050     if (fCutOnDThetaX && TMath::Abs(dThetaX)>fDThetaWindow) continue;
1051     //    double phi    = mlt->GetPhi(itr);
1052     double dphi   = mlt->GetDeltaPhi(itr);
1053     double dist   = mlt->CalcDist(itr);
1054     //
1055     FillHistosSet(histos,eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist);
1056     // special handling for mc info
1057     if (fillMC && fStack) {
1058       int lab0 = mlt->GetLabel(itr,0);
1059       int lab1 = mlt->GetLabel(itr,1);
1060       int typeMC = 2; // comb.bg.
1061       if (lab0 == lab1) typeMC = fStack->IsPhysicalPrimary(lab0) ? 0:1; // prim or sec
1062       if      (typeMC==0) FillHistosSet(fHistosTrPrim,eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist); // primary
1063       else if (typeMC==1) FillHistosSet(fHistosTrSec, eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist); // secondary
1064       else {
1065         FillHistosSet(fHistosTrComb,eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist); // comb
1066         // for combinatorals fill also the uncorrelated part
1067         if (fMultReco) {
1068           float *trl = fMultReco->GetTracklet(itr);
1069           int clId0 = (int)trl[AliITSMultReconstructor::kClID1];
1070           int clId1 = (int)trl[AliITSMultReconstructor::kClID2];
1071           float *clLabs0 = fMultReco->GetClusterOfLayer(0,clId0) + AliITSMultReconstructor::kClMC0;
1072           float *clLabs1 = fMultReco->GetClusterOfLayer(1,clId1) + AliITSMultReconstructor::kClMC0;
1073           if (!HaveCommonParent(clLabs0,clLabs1))
1074             FillHistosSet(fHistosTrCombU,eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist);
1075         }
1076       } // combinatorials
1077
1078       if (dist<fNStdCut) {
1079         double dphiS  = TMath::Abs(dphi) - fDPhiShift; if (dphi<0) dphiS = -dphiS;
1080         if (dphiS<fDPhiSCut) FillSpecies(typeMC, lab0);
1081       }
1082       if (fCheckReconstructables) CheckReconstructables();
1083     }
1084   }
1085   //
1086   //-------------------------------------------------------------TMP RS - singles ------->>>
1087   int offsH = fCurrCentBin*kNStandardH;
1088   TH2* hSingles = (TH2*)histos->UncheckedAt(offsH+kHEtaZvSPD1);
1089   if (hSingles) {
1090     int nclS = mlt->GetNumberOfSingleClusters();
1091     double *thtS = mlt->GetThetaSingle();
1092     for (int ics=nclS;ics--;) {
1093       double etaS = -TMath::Log(TMath::Tan(thtS[ics]/2));
1094       if (etaS<fEtaMin || etaS>fEtaMax) continue;
1095       hSingles->Fill(etaS,fESDVtx[2]);
1096     }
1097   }
1098   //-------------------------------------------------------------TMP RS - singles -------<<<
1099   //
1100 }
1101
1102 //_________________________________________________________________________
1103 void AliTrackletTaskMulti::FillMCPrimaries()
1104 {
1105   // fill all MC primaries Zv vs Eta
1106   if (!fStack || !fMCEvent) return;
1107
1108   //---------------------------------------- CHECK ------------------------------>>>
1109   TArrayF vtxMC;
1110   AliGenHijingEventHeader* pyHeader = 0;
1111   //
1112   if (fUseMC) {
1113     pyHeader = (AliGenHijingEventHeader*) fMCEvent->GenEventHeader();//header->GenEventHeader();
1114     pyHeader->PrimaryVertex(vtxMC);
1115   }
1116   //---------------------------------------- CHECK ------------------------------<<<
1117   //
1118   int ntr = fStack->GetNtrack();
1119   TH2* hprimEtaZ = (TH2F*)fHistosCustom->UncheckedAt(kHZVEtaPrimMC+fCurrCentBin);
1120   int nprim = 0;
1121   for (int itr=ntr;itr--;) {
1122     if (!fStack->IsPhysicalPrimary(itr)) continue;
1123     AliMCParticle *part  = (AliMCParticle*)fMCEvent->GetTrack(itr);
1124     if (!part->Charge()) continue;
1125     //
1126     //---------------------------------------- CHECK ------------------------------>>>
1127     /*
1128     Float_t dz = part->Zv() - vtxMC[2];
1129     if (TMath::Abs(dz)>1e-6) continue; // reject
1130     */
1131     //---------------------------------------- CHECK ------------------------------<<<
1132     //
1133     Float_t theta = part->Theta();
1134     if (theta<1e-6 || theta>TMath::Pi()-1e-6) continue;
1135     Float_t eta = part->Eta();
1136     if (eta<fEtaMin || eta>fEtaMax) continue;
1137     hprimEtaZ->Fill(eta, fESDVtx[2]);
1138     nprim++;
1139   }
1140   //
1141   ((TH1*)fHistosCustom->UncheckedAt(kHNPrimMeanMC))->Fill(fCurrCentBin,nprim);
1142   if (fNPart>0) {
1143     ((TH1*)fHistosCustom->UncheckedAt(kHNPrim2PartMC))->Fill(fCurrCentBin,nprim/fNPart);
1144     ((TH2*)fHistosCustom->UncheckedAt(kHNPrim2PartNpMC))->Fill(fNPart,nprim/fNPart);
1145     ((TH2*)fHistosCustom->UncheckedAt(kHNPartMC))->Fill(fNPart,fCurrCentBin);
1146     ((TH1*)fHistosCustom->UncheckedAt(kHNPartMeanMC))->Fill(fCurrCentBin,fNPart);
1147   }
1148   if (fNBColl>0) {
1149     ((TH1*)fHistosCustom->UncheckedAt(kHNPrim2BCollMC))->Fill(fCurrCentBin,nprim/fNBColl);
1150     ((TH2*)fHistosCustom->UncheckedAt(kHNPrim2BCollNpMC))->Fill(fNPart,nprim/fNBColl);
1151     ((TH2*)fHistosCustom->UncheckedAt(kHNBCollMC))->Fill(fNBColl,fCurrCentBin);
1152     ((TH1*)fHistosCustom->UncheckedAt(kHNBCollMeanMC))->Fill(fCurrCentBin,fNBColl);
1153   }
1154   //
1155 }
1156
1157 //_________________________________________________________________________
1158  void AliTrackletTaskMulti::FillHistosSet(TObjArray* histos, double eta,
1159                                           //double /*phi*/,double /*theta*/,
1160                                           double dphi,double dtheta,double dThetaX,
1161                                           double dist)
1162 {
1163   // fill standard set of histos
1164   if (dist>fNStdDev) return;
1165   //
1166   double dphiS  = TMath::Abs(dphi) - fDPhiShift;
1167   if (dphi<0) dphiS = -dphiS;
1168   //
1169   int offs = fCurrCentBin*kNStandardH;
1170   //
1171   if (histos->UncheckedAt(offs+kHDPhiSDThetaX))
1172     ((TH2*)histos->UncheckedAt(offs+kHDPhiSDThetaX))->Fill(dphiS,dThetaX);
1173   //
1174   if (histos->UncheckedAt(offs+kHDPhiDTheta))
1175     ((TH2*)histos->UncheckedAt(offs+kHDPhiDTheta))->Fill(dphi,dtheta);
1176   //
1177   if (histos->UncheckedAt(kHWDist))
1178     ((TH2*)histos->UncheckedAt(offs+kHWDist))->Fill(dist);
1179   //
1180   if (dist<fNStdCut && dphiS<fDPhiSCut && histos->UncheckedAt(offs+kHEtaZvCut))
1181     ((TH2*)histos->UncheckedAt(offs+kHEtaZvCut))->Fill(eta,fESDVtx[2]);
1182   //
1183 }
1184
1185 //__________________________________________________________________
1186 void AliTrackletTaskMulti::FillSpecies(Int_t primsec, Int_t id)
1187 {
1188   // fill PDGcode
1189   TH1 *hPart=0,*hParent=0;
1190   if (primsec==0) {
1191     hPart   = (TH1*)fHistosCustom->UncheckedAt(kHPrimPDG);
1192     hParent = (TH1*)fHistosCustom->UncheckedAt(kHPrimParPDG);
1193   }
1194   else if (primsec==1) {
1195     hPart   = (TH1*)fHistosCustom->UncheckedAt(kHSecPDG);
1196     hParent = (TH1*)fHistosCustom->UncheckedAt(kHSecParPDG);
1197   }
1198   else return;
1199   int ntr = fStack->GetNtrack();
1200   TParticle* part = fStack->Particle(id);
1201   int pdgCode = TMath::Abs(part->GetPdgCode());
1202   int pdgBin = GetPdgBin(pdgCode);
1203   int parID = part->GetFirstMother();
1204   int pdgCodePar = -1;
1205   int pdgBinPar = -1;
1206   while (parID>=0 && parID<ntr) {
1207     part = fStack->Particle(parID);
1208     pdgCodePar = TMath::Abs(part->GetPdgCode());
1209     parID = part->GetFirstMother();
1210   }
1211   if (pdgCodePar>0) pdgBinPar = GetPdgBin(pdgCodePar);
1212   //
1213   hPart->Fill(pdgBin,fCurrCentBin);
1214   hParent->Fill(pdgBinPar,fCurrCentBin);
1215   //
1216 }
1217
1218 //_________________________________________________________________________
1219 Int_t AliTrackletTaskMulti::GetCentralityBin(Float_t perc) const
1220 {
1221   // calculate centrality bin
1222   for (int i=0;i<fNCentBins;i++) if (perc>=fgkCentPerc[i] && perc<=fgkCentPerc[i+1]) return i;
1223   return -1;
1224 }
1225
1226 //_________________________________________________________________________
1227 Int_t AliTrackletTaskMulti::GetPdgBin(Int_t pdgCode)
1228 {
1229   // return my pdg bin
1230   int ncodes = sizeof(fgkPDGCodes)/sizeof(int);
1231   int pdgBin=0;
1232   for (pdgBin=0;pdgBin<ncodes;pdgBin++) if (pdgCode==fgkPDGCodes[pdgBin]) break;
1233   if (pdgBin>=ncodes) {
1234     if (float(pdgCode)>1e9) pdgBin = ncodes; // nuclei
1235     else pdgBin = ncodes+1; // unknown
1236   }
1237   return pdgBin;
1238 }
1239
1240 //_________________________________________________________________________
1241 Bool_t AliTrackletTaskMulti::HaveCommonParent(const float* clLabs0,const float* clLabs1)
1242 {
1243   // do 2 clusters have common parrent
1244   const int kMaxPar = 50;
1245   static int pars[2][50];
1246   int npars[2]={0,0};
1247   const float *labs[2] = {clLabs0,clLabs1};
1248   int ntr = fStack->GetNtrack();
1249   for (int il=0;il<2;il++) {
1250     for (int ilb=0;ilb<3;ilb++) {
1251       int lbl = (int)labs[il][ilb];
1252       if (lbl<0 || lbl>=ntr) continue;
1253       //
1254       while (npars[il]<kMaxPar-1) {
1255         pars[il][ npars[il]++ ] = lbl;
1256         TParticle* part = fStack->Particle(lbl);
1257         if (!part) break;
1258         lbl = part->GetFirstMother();
1259         if (lbl<1 || lbl>=ntr) break;
1260       }
1261     }
1262   }
1263   // compare array of parents
1264   for (int i0=npars[0];i0--;) for (int i1=npars[1];i1--;) if (pars[0][i0]==pars[1][i1]) return kTRUE;
1265   return kFALSE;
1266 }
1267
1268
1269 //_________________________________________________________________________
1270 void AliTrackletTaskMulti::CheckReconstructables()
1271 {
1272   // fill reconstructable tracklets hitsos
1273   static TArrayI trInd;
1274   static TBits   isPrimArr;
1275   //
1276   if (!fMultReco || !fMultReco->IsRecoDone()) {AliInfo("To check reconstructables the reco had to be requested"); return;}
1277   if (!fStack) {AliInfo("MC Stack is not availalble"); return;}
1278   const double kPtMin = 0.05;
1279   //
1280   TClonesArray *clArr[2];
1281   for (int ilr=0;ilr<2;ilr++) {
1282     clArr[ilr] = fMultReco->GetClustersOfLayer(ilr);
1283     if (!clArr[ilr]) {AliInfo("Clusters are not available"); return;}
1284   }
1285   //
1286   int ntr = fStack->GetNtrack();
1287   if (!ntr) return;
1288   trInd.Reset();
1289   if (trInd.GetSize()<ntr) trInd.Set(ntr);
1290   isPrimArr.ResetAllBits();
1291   // count track wich may be reconstructable
1292   //
1293   int ntrStore=0,ntrStorePrim=0;
1294   Int_t *trIndArr = trInd.GetArray();
1295   for (Int_t it=0; it<ntr; it++) {
1296     TParticle* part = fStack->Particle(it);
1297     if (TMath::Abs(part->Eta())>2.2)       continue;
1298     if (TMath::Abs(part->Pt())<kPtMin)      continue;
1299     if (fStack->IsPhysicalPrimary(it)) {
1300       isPrimArr.SetBitNumber(it);
1301       ntrStorePrim++;
1302     }
1303     else { // check if secondary is worth cheking
1304       TParticlePDG* pdgPart = part->GetPDG();
1305       if (TMath::Abs(pdgPart->Charge())!=3)  continue;
1306       if (part->R()>5.)                      continue;
1307     }
1308     trIndArr[it] = ++ntrStore;
1309   }
1310   //
1311   AliInfo(Form("Selected %d MC particles (%d prim) out of %d in the stack\n",ntrStore,ntrStorePrim,ntr));
1312   //
1313   const int kMaxCl=3;
1314   AliITSRecPoint **clIndL[2];
1315   clIndL[0] = new AliITSRecPoint*[kMaxCl*ntrStore]; // max 2 clusters per layer
1316   clIndL[1] = new AliITSRecPoint*[kMaxCl*ntrStore]; // max 2 clusters per layer
1317   memset(clIndL[0],0,kMaxCl*ntrStore*sizeof(AliITSRecPoint*));
1318   memset(clIndL[1],0,kMaxCl*ntrStore*sizeof(AliITSRecPoint*));
1319   //
1320   for (int ilr=0;ilr<2;ilr++) {
1321     TClonesArray *clusters = clArr[ilr];
1322     int ncl = clusters->GetEntriesFast();
1323     for (int icl=ncl;icl--;) {
1324       AliITSRecPoint *cl = (AliITSRecPoint*)clusters->UncheckedAt(icl);
1325       for (int ilb=3;ilb--;) {
1326         int lbl = cl->GetLabel(ilb); if (lbl<0 || lbl>=ntr) continue;
1327         int lblI = trIndArr[lbl];
1328         if (--lblI<0) continue;    // not kept
1329         for (int icc=0;icc<kMaxCl;icc++) if (!clIndL[ilr][lblI+icc*ntrStore]) {clIndL[ilr][lblI+ntrStore*icc] = cl; break;} // first empty one
1330       }
1331     }
1332   }
1333   //
1334   Float_t clusterLay[2][AliITSMultReconstructor::kClNPar];
1335   double trComp[6][kMaxCl*kMaxCl];
1336   int indQual[kMaxCl*kMaxCl];
1337   //
1338   for (int itr=ntr;itr--;) {
1339     int lblI = trIndArr[itr];
1340     if (--lblI<0) continue; // discarded
1341     int ntrCand = 0;        // number of tracklet candidates (including overlaps)
1342     for (int icl0=0;icl0<kMaxCl;icl0++) {
1343       AliITSRecPoint *cl0 = clIndL[0][lblI+icl0*ntrStore];
1344       if (!cl0 || !clIndL[1][lblI]) break;
1345       cl0->GetGlobalXYZ( clusterLay[0] );
1346       fMultReco->ClusterPos2Angles(clusterLay[0], fESDVtx);
1347       for (int icl1=0;icl1<kMaxCl;icl1++) {
1348         AliITSRecPoint *cl1 = clIndL[1][lblI+icl1*ntrStore];
1349         if (!cl1) break;
1350         cl1->GetGlobalXYZ( clusterLay[1] );
1351         fMultReco->ClusterPos2Angles(clusterLay[1], fESDVtx);
1352         trComp[AliITSMultReconstructor::kTrPhi][ntrCand]    = clusterLay[0][AliITSMultReconstructor::kClPh];
1353         trComp[AliITSMultReconstructor::kTrTheta][ntrCand]  = clusterLay[0][AliITSMultReconstructor::kClTh];
1354         trComp[AliITSMultReconstructor::kTrDTheta][ntrCand] = clusterLay[0][AliITSMultReconstructor::kClTh] - clusterLay[1][AliITSMultReconstructor::kClTh];
1355         trComp[AliITSMultReconstructor::kTrDPhi][ntrCand]   = clusterLay[0][AliITSMultReconstructor::kClPh] - clusterLay[1][AliITSMultReconstructor::kClPh];
1356         trComp[AliITSMultReconstructor::kTrLab1][ntrCand]   = icl1*10 + icl0;
1357         double &dphi = trComp[ntrCand][3];
1358         if (dphi>TMath::Pi()) dphi=2.*TMath::Pi()-dphi;     // take into account boundary condition
1359         trComp[5][ntrCand] = fMultReco->CalcDist(trComp[AliITSMultReconstructor::kTrDPhi][ntrCand],
1360                                                  trComp[AliITSMultReconstructor::kTrDTheta][ntrCand],
1361                                                  trComp[AliITSMultReconstructor::kTrTheta][ntrCand]);
1362         ntrCand++;
1363       }
1364     }
1365     if (!ntrCand) continue; // no tracklets
1366     if (ntrCand>1) TMath::Sort(ntrCand,trComp[5],indQual,kFALSE); else indQual[0] = 0; // sort in weighted distance
1367     if (fRemoveOverlaps) ntrCand = 1; // select the best
1368     //
1369     // disable worst tracklet with shared cluster
1370     for (int itc=0;itc<ntrCand;itc++) {
1371       int ind = indQual[itc];
1372       if (trComp[AliITSMultReconstructor::kTrLab1][ind]<0) continue; // already disabled
1373       for (int jtc=itc+1;jtc<ntrCand;jtc++) {
1374         int jnd = indQual[jtc];
1375         if (trComp[AliITSMultReconstructor::kTrLab1][jnd]<0) continue; // already disabled
1376         if ( int(trComp[AliITSMultReconstructor::kTrLab1][ind])/10 == int(trComp[AliITSMultReconstructor::kTrLab1][jnd])/10 ||
1377              int(trComp[AliITSMultReconstructor::kTrLab1][ind])%10 == int(trComp[AliITSMultReconstructor::kTrLab1][jnd])%10) trComp[AliITSMultReconstructor::kTrLab1][jnd] = -1;
1378       }
1379     }
1380     //
1381     // store, but forbid cluster reusing
1382     TObjArray* histos = isPrimArr.TestBitNumber(itr) ? fHistosTrRcblPrim : fHistosTrRcblSec;
1383     for (int itc=0;itc<ntrCand;itc++) {
1384       int ind = indQual[itc];
1385       if (trComp[4][ind]<0) continue; // discarded
1386       double eta    = -TMath::Log(TMath::Tan(trComp[AliITSMultReconstructor::kTrTheta][ind]/2));
1387       if (eta<fEtaMin || eta>fEtaMax) continue;
1388       double dThetaX = trComp[AliITSMultReconstructor::kTrTheta][ind];
1389       if (fScaleDTBySin2T) {
1390         double sint   =  TMath::Sin(trComp[AliITSMultReconstructor::kTrTheta][ind]);
1391         dThetaX /= (sint*sint);
1392       }
1393       FillHistosSet(histos,eta,
1394                     //trComp[AliITSMultReconstructor::kTrPhi][ind],trComp[AliITSMultReconstructor::kTrTheta][ind],
1395                     trComp[AliITSMultReconstructor::kTrDPhi][ind],trComp[AliITSMultReconstructor::kTrDTheta][ind],
1396                     dThetaX,trComp[5][ind]);
1397     }
1398   }
1399   //
1400   delete[] clIndL[0];
1401   delete[] clIndL[1];
1402 }
1403
1404 //_________________________________________________________________________
1405 void AliTrackletTaskMulti::FillClusterInfo()
1406 {
1407   // fill info on clusters associated to good tracklets
1408   if (!fMultReco) return;
1409   int ntr = fMultReco->GetNTracklets();
1410   int clID[2];
1411   TH2F *hclU[2] = {(TH2F*)fHistosCustom->UncheckedAt(kHClUsedInfoL0),(TH2F*)fHistosCustom->UncheckedAt(kHClUsedInfoL1)};
1412   TH2F *hclA[2] = {(TH2F*)fHistosCustom->UncheckedAt(kHClAllInfoL0),(TH2F*)fHistosCustom->UncheckedAt(kHClAllInfoL1)};
1413   for (int itr=ntr;itr--;) {
1414     Float_t *trc = fMultReco->GetTracklet(itr);
1415     if (TMath::Abs(TMath::Abs(trc[AliITSMultReconstructor::kTrDPhi])-fDPhiShift)>fDPhiSCut) continue;
1416     if (fMultReco->CalcDist(trc[AliITSMultReconstructor::kTrDPhi],
1417                             trc[AliITSMultReconstructor::kTrDTheta],
1418                             trc[AliITSMultReconstructor::kTrTheta]) > fNStdCut) continue;
1419     clID[0] = (int)trc[AliITSMultReconstructor::kClID1];
1420     clID[1] = (int)trc[AliITSMultReconstructor::kClID2];
1421     for (int il=0;il<2;il++) {
1422       Float_t *clinf = fMultReco->GetClusterOfLayer(il,clID[il]);
1423       hclU[il]->Fill( clinf[AliITSMultReconstructor::kClZ], clinf[AliITSMultReconstructor::kClPh]);
1424     }
1425   }
1426   //
1427   for (int il=0;il<2;il++) for (int ic=fMultReco->GetNClustersLayer(il);ic--;) {
1428       Float_t *clinf = fMultReco->GetClusterOfLayer(il,ic);
1429       hclA[il]->Fill( clinf[AliITSMultReconstructor::kClZ], clinf[AliITSMultReconstructor::kClPh]);
1430     }
1431   //
1432 }
1433
1434 //_________________________________________________________________________
1435 void AliTrackletTaskMulti::FillClusterInfoFromMult(const AliMultiplicity* mlt, double zVertex)
1436 {
1437   // fill info on clusters taking them from Multiplicity object
1438   const double kRSPD2 = 3.9;
1439   TH2F *hclU = (TH2F*)fHistosCustom->UncheckedAt(kHClUsedInfoL0);
1440   TH2F *hclA = (TH2F*)fHistosCustom->UncheckedAt(kHClAllInfoL0);
1441   int ntr = mlt->GetNumberOfTracklets();
1442   for (int itr=ntr;itr--;) {
1443     Bool_t goodTracklet = kTRUE;
1444     if (TMath::Abs( mlt->GetDeltaPhi(itr)-fDPhiShift)>fDPhiSCut) goodTracklet = kFALSE;
1445     if (mlt->CalcDist(itr) > fNStdCut) goodTracklet = kFALSE;
1446     double phi   = mlt->GetPhi(itr);
1447     double z     = kRSPD2/TMath::Tan(mlt->GetTheta(itr)) + zVertex;
1448     if (goodTracklet) hclU->Fill(z,phi);
1449     hclA->Fill(z,phi);
1450   }
1451   //
1452   int ncl = mlt->GetNumberOfSingleClusters();
1453   for (int icl=ncl;icl--;) {
1454     double phi   = mlt->GetPhiSingle(icl);
1455     double z     = kRSPD2/TMath::Tan(mlt->GetThetaSingle(icl)) + zVertex;
1456     hclA->Fill(z,phi);
1457   }
1458   //
1459 }
1460