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