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