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