AliPhysics  58f3d52 (58f3d52)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliTrackletTaskMulti.cxx
Go to the documentation of this file.
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 
17 // Class AliTrackletTaskMulti //
18 // Analysis task to produce data and MC histos needed for tracklets //
19 // dNdEta extraction in multiple bins in one go //
20 // Author: ruben.shahoyan@cern.ch //
21 // Hacker: roberto.preghenella@bo.infn.it (R+) //
23 /*
24  Important parameters to set:
25  1) make sure to initialize correct geometry in UserCreateOutputObjects
26  2) The cut on signal selection variable (delta, dphi ...) should be decided beforehand
27 ...
28 */
29 
30 #include "TFile.h"
31 #include "TChain.h"
32 #include "TTree.h"
33 #include "TRandom.h"
34 #include "TH1F.h"
35 #include "TH2F.h"
36 #include "TH3F.h"
37 #include "THnSparse.h"
38 #include "TList.h"
39 #include "TNtuple.h"
40 #include "TObjArray.h"
41 #include "TGeoGlobalMagField.h"
42 
43 #include "AliAnalysisManager.h"
44 
45 #include "AliMultiplicity.h"
46 #include "AliESDEvent.h"
47 #include "AliESDInputHandler.h"
48 #include "AliESDInputHandlerRP.h"
49 #include "AliCDBPath.h"
50 #include "AliCDBManager.h"
51 #include "AliCDBEntry.h"
52 #include "AliCDBStorage.h"
53 #include "AliGeomManager.h"
54 #include "AliMagF.h"
55 #include "AliESDVZERO.h"
56 #include "AliESDZDC.h"
57 #include "AliRunLoader.h"
58 #include "AliMCEventHandler.h"
59 #include "AliMCEvent.h"
60 #include "AliMCParticle.h"
61 #include "AliStack.h"
62 #include "AliGenEventHeader.h"
63 #include "AliCentrality.h"
64 #include "AliTriggerAnalysis.h"
65 #include "AliMultSelection.h"
66 #include "AliITSsegmentationSPD.h"
67 #include "AliITSRecPoint.h"
68 #include "AliITSgeomTGeo.h"
69 #include "AliITSMultReconstructor.h"
70 
71 #include "AliLog.h"
72 
73 #include "AliPhysicsSelection.h"
74 #include "AliTrackletTaskMulti.h"
75 #include "AliITSMultRecBg.h"
76 #include "AliGenEventHeader.h"
77 #include "AliGenHijingEventHeader.h"
78 #include "AliGenDPMjetEventHeader.h"
79 #include "AliGenCocktailEventHeader.h"
80 #include "AliESDtrackCuts.h"
81 
82 ClassImp(AliTrackletTaskMulti)
83 
84 
85 // 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14
86 const char* AliTrackletTaskMulti::fgCentSelName[] = {"V0M","V0A","V0C","FMD","TRK","TKL","CL0","CL1","V0MvsFMD","ZNA","TKLvsV0M","ZEMvsZDC","V0A123","V0A0","V0S","MB","V0Mplus05","V0Mminus05","SPDClustersCorr"};
87 
88 const char* AliTrackletTaskMulti::fgkPDGNames[] = {
89 "#pi^{+}",
90 "p",
91 "K^{+}",
92 "K^{*+}",
93 "e^{-}",
94 "#mu^{-}",
95 "#rho^{+}",
96 "D^{+}",
97 "D^{*+}",
98 "D_{s}^{+}",
99 "D_{s}^{*+}",
100 "#Delta^{-}",
101 "#Delta^{+}",
102 "#Delta^{++}",
103 "#Sigma^{-}",
104 "#Sigma^{+}",
105 "#Sigma^{*-}",
106 "#Sigma^{*+}",
107 "#Sigma^{*+}_{c}",
108 "#Sigma^{*++}_{c}",
109 "#Xi^{-}",
110 "#Xi^{*-}",
111 "#Lambda^{+}_{c}",
112 "n",
113 "#Delta^{0}",
114 "#gamma",
115 "K^{0}_{S}",
116 "K^{0}_{L}",
117 "K^{0}",
118 "K^{*}",
119 "#eta",
120 "#pi^{0}",
121 "#rho^{0}",
122 "#varphi",
123 "#eta'",
124 "#omega",
125 "#Lambda",
126 "#Sigma^{0}",
127 "#Sigma^{*0}_{c}",
128 "#Sigma^{*0}",
129 "D^{0}",
130 "D^{*0}",
131 "#Xi_{0}",
132 "#Xi^{*0}",
133 "#Xi^{0}_{c}",
134 "#Xi^{*0}_{c}",
135 "Nuclei",
136 "Others"
137 };
138 
139 const int AliTrackletTaskMulti::fgkPDGCodes[] = {
140  211,
141  2212,
142  321,
143  323,
144  11,
145  13,
146  213,
147  411,
148  413,
149  431,
150  433,
151  1114,
152  2214,
153  2224,
154  3112,
155  3222,
156  3114,
157  3224,
158  4214,
159  4224,
160  3312,
161  3314,
162  4122,
163  2112,
164  2114,
165  22,
166  310,
167  130,
168  311,
169  313,
170  221,
171  111,
172  113,
173  333,
174  331,
175  223,
176  3122,
177  3212,
178  4114,
179  3214,
180  421,
181  423,
182  3322,
183  3324,
184  4132,
185  4314
186 // nuclei
187 // unknown
188 };
189 
190 const double kEtaBinWidth=0.25;
191 
192 //________________________________________________________________________
193 /*//Default constructor
194 AliTrackletTaskMulti::AliTrackletTaskMulti(const char *name)
195  : AliAnalysisTaskSE(name),
196 */
197 //________________________________________________________________________
199  : AliAnalysisTaskSE(name),
200  //
201  fOutput(0),
202  //
203  fDoNormalReco(kFALSE),
204  fDoInjection(kFALSE),
205  fDoRotation(kFALSE),
206  //
207  fUseMC(kFALSE),
208  fCheckReconstructables(kFALSE),
209  //
210  fHistosTrData(0),
211  fHistosTrInj(0),
212  fHistosTrRot(0),
213  //
214  fHistosTrPrim(0),
215  fHistosTrSec(0),
216  fHistosTrComb(0),
217  fHistosTrCombU(0),
218  //
219  fHistosTrRcblPrim(0),
220  fHistosTrRcblSec(0),
221  fHistosCustom(0),
222  //
223  fEtaMin(-3.0),
224  fEtaMax(3.0),
225  fPhiMin(-999.),
226  fPhiMax(999.),
227  fZVertexMin(-20),
228  fZVertexMax( 20),
229  //
230  fScaleDTBySin2T(kFALSE),
231  fCutOnDThetaX(kFALSE),
232  fNStdDev(1.),
233  fDPhiWindow(0.08),
234  fDThetaWindow(0.025),
235  fDPhiShift(0.0045),
236  fPhiOverlapCut(0.005),
237  fZetaOverlap(0.05),
238  fPhiRot(0.),
239  fInjScale(1.),
240  fRemoveOverlaps(kFALSE),
241  //
242  fDPhiSCut(0.06),
243  fNStdCut(1.),
244  fMCV0Scale(1.),
245  //
246  fTrigSel(AliVEvent::kINT7),
247  //
248  fMultReco(0),
249  fRPTree(0),
250  fStack(0),
251  fMCEvent(0),
252  //
253  fNPrimMCeta2(0),
254  fNTreta2(0),
255  fNPart(0),
256  fNBColl(0),
257  fCurrCentBin(-1),
258  fNCentBins(0),
259  fCentPerc(),
260  fUseCentralityVar("V0M"),
261  fIsSelected(kFALSE),
262  fVtxOK(kFALSE),
263  fUseSpecialOutput(kFALSE),
264  fReweightStack(0x0), // R+HACK
265  fReweightFlag(0), // R+HACK
266  fEtaBinWidth(0.25) // cholm
267 {
268  // Constructor
269 
270  DefineOutput(1, TList::Class());
271  //
273  SetNStdDev();
274  SetPhiWindow();
275  SetThetaWindow();
276  SetPhiShift();
279  SetPhiRot();
281  //
282  fCentPerc.Set(2);
283  fCentPerc[0] = 0.0;
284  fCentPerc[1] = 100.0;
285  //
286 }
287 
288 //________________________________________________________________________
290 {
291  // Destructor
292  // histograms are in the output list and deleted when the output
293  // list is deleted by the TSelector dtor
294  if (fOutput && !AliAnalysisManager::GetAnalysisManager()->IsProofMode()) { //RRR
295  printf("Deleteing output\n");
296  delete fOutput;
297  fOutput = 0;
298  }
299  //
300  delete fMultReco;
301  //
302  delete fHistosTrData;
303  delete fHistosTrPrim;
304  delete fHistosTrSec;
305  delete fHistosTrComb;
306  delete fHistosTrCombU;
307  delete fHistosTrInj;
308  delete fHistosTrRot;
309  delete fHistosTrRcblPrim;
310  delete fHistosTrRcblSec;
311  delete fHistosCustom;
312  //
313 }
314 
315 //________________________________________________________________________
317 {
318  Printf("%s - %s", GetName(), GetTitle());
319  TH1* hstat;
320  TList* lst = dynamic_cast<TList*>(GetOutputData(1));
321  if (!lst) AliWarning("No output list yet");
322  else {
323  TH1* hstat=(TH1*)lst->FindObject("hStat");
324  if (!hstat)
325  AliWarning("No statistics histogram yet");
326  else {
327  TH1* stats = static_cast<TH1*>(hstat->Clone());
328  stats->SetDirectory(0);
329  TAxis* axis = stats->GetXaxis();
330  for (Int_t i = 1; i <= 4; i++) {
331  Printf("%31s: %f", axis->GetBinLabel(i), stats->GetBinContent(i));
332  }
333  Double_t scale = stats->GetBinContent(4);
334  stats->Scale(1/scale);
335  for (Int_t i = 6; i <= 23; i++) {
336  Printf("%31s: %f", axis->GetBinLabel(i), stats->GetBinContent(i));
337  }
338  delete stats;
339  }
340  }
341  Printf("--- Overall options ---");
342  Printf("DoNormalReco: %s", (fDoNormalReco ? "yes" : "no"));
343  Printf("DoInjection: %s", (fDoInjection ? "yes" : "no"));
344  Printf("DoRotation: %s", (fDoRotation ? "yes" : "no"));
345  Printf("UseMC: %s", (fUseMC ? "yes" : "no"));
346  Printf("CheckReconstructables: %s", (fCheckReconstructables ?
347  "yes" : "no"));
348  Printf("--- Tracklet reco settings ---");
349  Printf("EtaMin: %f", fEtaMin);
350  Printf("EtaMax: %f", fEtaMax);
351  Printf("PhiMin: %f", fPhiMin);
352  Printf("PhiMax: %f", fPhiMax);
353  Printf("ZVertexMin: %f", fZVertexMin);
354  Printf("ZVertexMax: %f", fZVertexMax);
355  Printf("--- Tracklet cuts --- ");
356  Printf("ScaleDTBySin2T: %s", (fScaleDTBySin2T ? "yes" : "no"));
357  Printf("CutOnDThetaX: %s", (fCutOnDThetaX ? "yes" : "no"));
358  Printf("NStdDev: %f", fNStdDev);
359  Printf("DPhiWindow: %f", fDPhiWindow);
360  Printf("DThetaWindow: %f", fDThetaWindow);
361  Printf("DPhiShift: %f", fDPhiShift);
362  Printf("PhiOverlapCut: %f", fPhiOverlapCut);
363  Printf("ZetaOverlap: %f", fZetaOverlap);
364  Printf("PhiRot: %f", fPhiRot);
365  Printf("InjScale: %f", fInjScale);
366  Printf("RemoveOverlaps: %s", (fRemoveOverlaps ? "yes" : "no"));
367  Printf("--- Task options ---");
368  Printf("DPhiSCut: %f", fDPhiSCut);
369  Printf("NStdCut: %f", fNStdCut);
370  Printf("MCV0Scale: %f", fMCV0Scale);
371  Printf("TrigSel: 0x%x", fTrigSel);
372  Printf("UseCentralityVar: %s", fUseCentralityVar.Data());
373  Printf("UseSpecialOutput: %s", (fUseSpecialOutput ? "yes" : "no"));
374  Printf("--- Reweight of stack ---");
375  Printf("ReweightStack: %d", fReweightStack);
376  Printf("ReweightFlag: %d", fReweightFlag);
377  Printf("EtaBinWidth: %f", fEtaBinWidth);
378 
379 }
380 
381 //________________________________________________________________________
383 {
384  //
385  if (fUseSpecialOutput) OpenFile(1);
386  fOutput = new TList();
387  fOutput->SetOwner();
388  //
389  //
390  Bool_t needGeom = GetDoNormalReco() || GetDoInjection() || GetDoRotation();
391  if (needGeom) {
392  AliCDBManager *man = AliCDBManager::Instance();
393  if (fUseMC) {
394  man->SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Residual");
395  AliCDBEntry* obj = man->Get("GRP/Geometry/Data",137161);
396  AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
397  if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",137161,-1,-1)) AliFatal("Failed to misalign geometry");
398  }
399  else {
400  man->SetDefaultStorage("alien://Folder=/alice/data/2010/OCDB"); //man->SetRun(137366);
401  AliCDBEntry* obj = man->Get("GRP/Geometry/Data",137161);
402  AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
403 
404  /* R+TEMP
405  man->UnsetDefaultStorage();
406  man->SetDefaultStorage("alien://Folder=/alice/cern.ch/user/s/shahoian/ITS_fixSPD72");
407  */
408 
409  if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",137161,-1,-1)) AliFatal("Failed to misalign geometry");
410 
411  }
412  }
413  //
414  // Create histograms
415  //---------------------------------------------Standard histos per tracklet type--->>
416  UInt_t hPattern = 0xffffffff;
417  fHistosTrData = BookHistosSet("TrData",hPattern);
418  // hPattern &= ~(BIT(kHEtaZvSPD1)); // fill single clusters for "data" only
419  if (GetDoInjection()) fHistosTrInj = BookHistosSet("TrInj",hPattern);
420  if (GetDoRotation()) fHistosTrRot = BookHistosSet("TrRot",hPattern);
421  if (fUseMC) {
422  fHistosTrPrim = BookHistosSet("TrPrim",hPattern);
423  fHistosTrSec = BookHistosSet("TrSec",hPattern);
424  fHistosTrComb = BookHistosSet("TrComb",hPattern);
425  fHistosTrCombU = BookHistosSet("TrCombU",hPattern);
427  fHistosTrRcblPrim = BookHistosSet("TrRcblPrim",hPattern);
428  fHistosTrRcblSec = BookHistosSet("TrRcblSec",hPattern);
429  }
430  }
431  //---------------------------------------------Standard histos per tracklet type---<<
432  //
433  //---------------------------------------------Custom Histos----------------------->>
434  // put here any non standard histos
436  //
437  //---------------------------------------------Custom Histos-----------------------<<
438  int nhist = fOutput->GetEntries();
439  for (int i=0;i<nhist;i++) {
440  TObject* hst = fOutput->At(i);
441  if (!hst || !(hst->InheritsFrom(TH1::Class()))) continue;
442  ((TH1*)hst)->Sumw2();
443  }
444  //
445  AliInfo(Form("Centrality type selected: %s\n",fUseCentralityVar.Data()));
446  PostData(1, fOutput);
447  //
448 }
449 
450 //________________________________________________________________________
452 {
453  // Main loop
454  //
455  static Bool_t statOK = kFALSE;
456  if (!statOK) {
457  RegisterStat();
458  statOK = kTRUE;
459  }
460  //
461  printf("Do: %d %d %d\n",GetDoNormalReco(),GetDoInjection(),GetDoRotation());
462  Bool_t needRecPoints = GetDoNormalReco() || GetDoInjection() || GetDoRotation();
463  //
464  AliAnalysisManager* anMan = AliAnalysisManager::GetAnalysisManager();
465  fRPTree = 0;
466  AliESDInputHandler *handler = (AliESDInputHandler*)anMan->GetInputEventHandler();
467  AliESDInputHandlerRP *handRP = 0;
468  if (needRecPoints &&
469  !handler->IsA()->InheritsFrom(AliESDInputHandlerRP::Class())) {
470  printf("No RP handler\n");
471  return;
472  }
473  if (needRecPoints) {
474  handRP = static_cast<AliESDInputHandlerRP*>(handler);
475  if (!handRP) { printf("No RP handler\n"); return; }
476  }
477  AliESDEvent *esd = static_cast<AliESDEvent*>(handler->GetEvent());
478  if (!esd) { printf("No AliESDEvent\n"); return; }
479  //
480  //
481  if (needRecPoints) {
482  Printf("Getting rec-points (clusters) from handler (%p)", handRP);
483  fRPTree = handRP->GetTreeR("ITS");
484  if (!fRPTree) { printf("Invalid ITS cluster tree !\n"); return; }
485  }
486  //
487  // do we need to initialize the field?
488  AliMagF* field = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
489  if (!field && !esd->InitMagneticField()) {printf("Failed to initialize the B field\n");return;}
490  //
491  AliMCEventHandler* eventHandler = 0;
492  fMCEvent = 0;
493  fStack = 0;
494  //
495  // ========================== MC EVENT INFO ====================================>>>
496  if (fUseMC) {
497  eventHandler = (AliMCEventHandler*)anMan->GetMCtruthEventHandler();
498  if (!eventHandler) { printf("ERROR: Could not retrieve MC event handler\n"); return; }
499  fMCEvent = eventHandler->MCEvent();
500  if (!fMCEvent) { printf("ERROR: Could not retrieve MC event\n"); return; }
501  fStack = fMCEvent->Stack();
502  if (!fStack) { printf("Stack not available\n"); return; }
503  }
504  //
505  // =============================================================================>>>
506  //
507  TH1* hstat = (TH1*)fHistosCustom->UncheckedAt(kHStat);
508  hstat->Fill(kEvTot0); // RS
509  //
510  // MC Generator info
511  AliGenEventHeader* mcGenH = 0;
512  AliGenHijingEventHeader* hHijing=0;
513  AliGenDPMjetEventHeader* hDpmJet=0;
514  fNPart = 0;
515  fNBColl = 0;
516  if (fUseMC) {
517  mcGenH = fMCEvent->GenEventHeader();
518  //
519  if (mcGenH->InheritsFrom(AliGenCocktailEventHeader::Class())) {
520  TList* headers = ((AliGenCocktailEventHeader*)mcGenH)->GetHeaders();
521  TIter next(headers);
522  AliGenEventHeader* mcGenH1 = 0;
523  while ( (mcGenH1=(AliGenEventHeader*)next()) ) {
524  if (mcGenH1->InheritsFrom(AliGenHijingEventHeader::Class())) {
525  hHijing = (AliGenHijingEventHeader*)mcGenH1;
526  break;
527  }
528  if (mcGenH1->InheritsFrom(AliGenDPMjetEventHeader::Class())) {
529  hDpmJet = (AliGenDPMjetEventHeader*)mcGenH1;
530  break;
531  }
532  }
533  }
534  //
535  if ( hHijing || mcGenH->InheritsFrom(AliGenHijingEventHeader::Class())) {
536  if (!hHijing) hHijing = (AliGenHijingEventHeader*)mcGenH;
537  fNPart = (hHijing->ProjectileParticipants()+hHijing->TargetParticipants())/2.;
538  fNBColl = hHijing->NN()+hHijing->NNw()+hHijing->NwN()+hHijing->NwNw();
539  }
540  else if ( hDpmJet || mcGenH->InheritsFrom(AliGenDPMjetEventHeader::Class())) {
541  if (!hDpmJet) hDpmJet = (AliGenDPMjetEventHeader*)mcGenH;
542  fNPart = (hDpmJet->ProjectileParticipants()+hDpmJet->TargetParticipants())/2.;
543  fNBColl = hDpmJet->NN()+hDpmJet->NNw()+hDpmJet->NwN()+hDpmJet->NwNw();
544  // Special for SD rejection
545  int npProj = hDpmJet->ProjectileParticipants(); // A
546  // int npTarg = hDpmJet->TargetParticipants(); // p
547  int nsd1,nsd2,ndd;
548  hDpmJet->GetNDiffractive(nsd1,nsd2,ndd);
549  if (ndd==0 && (npProj==nsd1 || npProj==nsd2)) {
550  hstat->Fill(kNEvSDMC);
551  return; // reject SD
552  }
553  }
554  else {} // unknown generator
555  //
556  TArrayF vtmc(3);
557  mcGenH->PrimaryVertex(vtmc);
558  for (int i=3;i--;) fVtxMC[i] = vtmc[i];
559  }
560  // =============================================================================>>>
561  //
562  Int_t trg = handler->IsEventSelected();
563  fIsSelected = trg & fTrigSel; //AliVEvent::kINT7;
564  if (fTrigSel == AliVEvent::kAny)
565  fIsSelected = kTRUE;
566  /*
567  if (!fUseMC) {
568  TString trigStr(esd->GetFiredTriggerClasses());
569  if (!trigStr.Contains("CINT5-B-")) return;
570  }
571  AliVVZERO* esdVZERO = esd->GetVZEROData();
572  fIsSelected = !((esdVZERO->GetV0ADecision()!=1) || (esdVZERO->GetV0CDecision()!=1));
573  */
574  /*
575  fIsSelected = kTRUE;
576  AliVVZERO* esdVZERO = esd->GetVZEROData();
577  if ((esdVZERO->GetV0ADecision()!=1) || (esdVZERO->GetV0CDecision()==2)) fIsSelected = kFALSE;
578  if (!fUseMC) {
579  TString trigStr(esd->GetFiredTriggerClasses());
580  if (!trigStr.Contains("CINT5-B-")) fIsSelected = kFALSE;
581  }
582  */
583  //
584  const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
585  const AliMultiplicity* multESD = esd->GetMultiplicity();
586  //
587  /* R+COMMENT
588  * OLD CENTRALITY SELECTION
589 
590  AliCentrality *centrality = esd->GetCentrality();
591  Double_t centPercentile = centrality->GetCentralityPercentileUnchecked(fUseCentralityVar.Data());
592  if (fUseMC && centPercentile==0&&multESD->GetNumberOfITSClusters(0)<1 &&
593  multESD->GetNumberOfITSClusters(1)<1) {
594  centPercentile = 100;
595  printf("Overriding centrality 0 to 100:\n"
596  "Cent%s:%.2f Sel:%d Ntrk:%d SPD1:%d SPD2:%d V0A:%f V0C:%f CentV0A:%.2f CentV0M:%.2f CentCl1:%.2f\n",
597  fUseCentralityVar.Data(),centPercentile,fIsSelected,multESD->GetNumberOfTracklets(),
598  multESD->GetNumberOfITSClusters(0),multESD->GetNumberOfITSClusters(1),
599  esd->GetVZEROData()->GetMTotV0A(),esd->GetVZEROData()->GetMTotV0C(),
600  centrality->GetCentralityPercentileUnchecked("V0A"),
601  centrality->GetCentralityPercentileUnchecked("V0M"),
602  centrality->GetCentralityPercentileUnchecked("CL1"));
603  }
604  */
605 
606  AliMultSelection *multSelection = (AliMultSelection *)esd->FindListObject("MultSelection");
607  Double_t centPercentile = multSelection->GetMultiplicityPercentile(fUseCentralityVar.Data());
608  //
609  if (fUseCentralityVar.EqualTo("MB"))
610  centPercentile = 0.;
611 
612  ((TH1*)fHistosCustom->UncheckedAt(kHCentDistNoSel))->Fill(centPercentile);
613  //
614  /*
615  if (centPercentile<1 || fIsSelected) {
616  printf("Cent%s:%.2f Sel:%d Ntrk:%d SPD1:%d SPD2:%d V0A:%f V0C:%f CentV0A:%.2f CentV0M:%.2f CentCl1:%.2f\n",
617  fUseCentralityVar.Data(),centPercentile,fIsSelected,multESD->GetNumberOfTracklets(),
618  multESD->GetNumberOfITSClusters(0),multESD->GetNumberOfITSClusters(1),
619  esd->GetVZEROData()->GetMTotV0A(),esd->GetVZEROData()->GetMTotV0C(),
620  centrality->GetCentralityPercentileUnchecked("V0A"),
621  centrality->GetCentralityPercentileUnchecked("V0M"),
622  centrality->GetCentralityPercentileUnchecked("CL1"));
623  }
624  */
625  const double kSafeMargin = 1e-3;
626  if (centPercentile<-kSafeMargin) return;
627 
628  if (centPercentile<kSafeMargin) centPercentile = kSafeMargin;
629  if (centPercentile>100-kSafeMargin) centPercentile = 100.-kSafeMargin;
630  //
631  fCurrCentBin = GetCentralityBin(centPercentile);
632  // R+REWEIGHT
633  if (fUseMC && fReweightStack) { printf(">>> Reweighting stack\n"); ReweightStack(1.); } // R+HACK
634  //
635  // ==================== STORE SOME GLOBAL INFO FOR ALL EVENTS ==============>>>
636  Float_t multV0=0;
637  AliESDVZERO* esdV0 = esd->GetVZEROData();
638  if (esdV0) {
639  multV0 = esdV0->GetMTotV0A()+esdV0->GetMTotV0C();
640  if (fUseMC) multV0 *= fMCV0Scale;
641  }
642  ((TH1*)fHistosCustom->UncheckedAt(kHV0NoSel))->Fill(multV0);
643  //
644  float nSPD2 = multESD->GetNumberOfITSClusters(1);
645  ((TH1*)fHistosCustom->UncheckedAt(kHNClSPD2NoSel))->Fill(nSPD2);
646  //
647  //------------------------------------------------------
648  AliESDZDC *esdZDC = esd->GetESDZDC();
649  float zdcEnergy=0,zemEnergy=0;
650  if (esdZDC) {
651  zdcEnergy = (esdZDC->GetZDCN1Energy() + esdZDC->GetZDCP1Energy() + esdZDC->GetZDCN2Energy()+ esdZDC->GetZDCP2Energy());
652  zemEnergy = (esdZDC->GetZDCEMEnergy(0)+ esdZDC->GetZDCEMEnergy(1))/8.;
653  }
654  // PutZDCSelection
655  ((TH2*)fHistosCustom->UncheckedAt(kHZDCZEMNoSel))->Fill(zemEnergy,zdcEnergy);
656  //
657  // ==================== STORE SOME GLOBAL INFO FOR ALL EVENTS ==============<<<
658  //
659  // printf("CentPerc: %f : Bin %d ZDC: %f ZEM: %f\n",centPercentile, fCurrCentBin, zdcEnergy,zemEnergy);
660  if (fCurrCentBin<0) {
661  //printf("Reject: %.1f : V0:%.1f V0Cor:%.1f V0CR:%.1f SPD2c:%.1f\n",mltTst, multV0,multV0Corr,multV0CorrResc,nSPD2Corr);
662  return;
663  }
665  //
666  fVtxOK = kFALSE;
667  for (int i=3;i--;) fESDVtx[i] = 0;
668  if (vtxESD->GetNContributors()>0) {
669  TString vtxTyp = vtxESD->GetTitle();
670  if ( !vtxTyp.Contains("vertexer: Z") || (vtxESD->GetDispersion()<0.04 && vtxESD->GetZRes()<0.25)) {
671  fVtxOK = kTRUE;
672  fESDVtx[0] = vtxESD->GetX();
673  fESDVtx[1] = vtxESD->GetY();
674  fESDVtx[2] = vtxESD->GetZ();
675  }
676  }
677  //
678  if (fIsSelected) hstat->Fill(kBinEntries+kEvPassPS + kEntriesPerBin*fCurrCentBin);
679  //
680  if (fVtxOK && fIsSelected) {
681  ((TH1*)fHistosCustom->UncheckedAt(kHZVtxNoSel))->Fill(fESDVtx[2]);
682  hstat->Fill(kBinEntries+kEvPassVtx + kEntriesPerBin*fCurrCentBin);
683  }
684  //
685  fVtxOK &= (fESDVtx[2] >= fZVertexMin && fESDVtx[2] <= fZVertexMax);
686  //
687  // if (!fVtxOK || !fIsSelected) return;
688 
689  if (fUseMC) {
690  FillMCPrimaries();
691  ((TH2*)fHistosCustom->UncheckedAt(kHZVtxMCNoPhSel))->Fill(fVtxMC[2],fCurrCentBin);
692  if (fIsSelected) ((TH2*)fHistosCustom->UncheckedAt(kHZVtxMCNoVtSel))->Fill(fVtxMC[2],fCurrCentBin);
693  if ((fVtxMC[2] < fZVertexMin || fVtxMC[2] > fZVertexMax)) hstat->Fill(kBinEntries+kEvInMltBin + kEntriesPerBin*fCurrCentBin);
694  //
695  }
696  if (!fVtxOK || !fIsSelected || centPercentile > 80.) return; // R+HACK
697  //
698  // ===================== Store multiplicity estimators ===============================
699  double etaRange = TMath::Abs(fEtaMax-fEtaMin)/2.;
700  double mltE0 = AliESDtrackCuts::GetReferenceMultiplicity(esd, AliESDtrackCuts::kTrackletsITSTPC, etaRange);
701  double mltE1 = AliESDtrackCuts::GetReferenceMultiplicity(esd, AliESDtrackCuts::kTrackletsITSSA , etaRange);
702  double mltE2 = AliESDtrackCuts::GetReferenceMultiplicity(esd, AliESDtrackCuts::kTracklets, etaRange);
703  ((TH1F*)fHistosCustom->UncheckedAt(kHMltEstTrITSTPC))->Fill(mltE0);
704  ((TH1F*)fHistosCustom->UncheckedAt(kHMltEstTrITSSA))->Fill(mltE1);
705  ((TH1F*)fHistosCustom->UncheckedAt(kHMltEstTr))->Fill(mltE2);
706  //
707  // ===================== Process event passing all selections ===============================
708  //
709  ((TH1*)fHistosCustom->UncheckedAt(kHCentDist))->Fill(centPercentile);
710  hstat->Fill(kEvTot); // RS
711  ((TH1*)fHistosCustom->UncheckedAt(kHStatCentBin))->Fill(fCurrCentBin);
712  ((TH1*)fHistosCustom->UncheckedAt(kHStatCent))->Fill(centPercentile);
713  //
714  // register Ntracklets and ZVertex of the event
715  if (fUseMC) ((TH2*)fHistosCustom->UncheckedAt(kHZVtxMC))->Fill(fVtxMC[2],fCurrCentBin);
716  ((TH2*)fHistosCustom->UncheckedAt(kHZVtx))->Fill(fESDVtx[2],fCurrCentBin);
717  ((TH2*)fHistosCustom->UncheckedAt(kHV0))->Fill(multV0,fCurrCentBin);
718  ((TH2*)fHistosCustom->UncheckedAt(kHNClSPD2))->Fill(nSPD2,fCurrCentBin);
719  ((TH3*)fHistosCustom->UncheckedAt(kHZDCZEM))->Fill(zemEnergy,zdcEnergy,fCurrCentBin);
720  //
721  // normal reconstruction
722  hstat->Fill(kBinEntries+kEvProcData + kEntriesPerBin*fCurrCentBin);
723  //
724  if (GetDoNormalReco() || GetDoInjection()) { // for the injection the normal reco should be done
725  InitMultReco();
726  fMultReco->Run(fRPTree, fESDVtx);
727  printf("Multiplicity Reconstructed:\n");
728  AliMultiplicity* mlt = fMultReco->GetMultiplicity();
729  if (mlt) mlt->Print();
730  if (GetDoNormalReco()) FillHistos(kData,mlt);
731  FillClusterInfo();
732  //
733  }
734  if (!GetDoNormalReco()) {
735  FillHistos(kData,multESD); // fill data histos from ESD
736  FillClusterInfoFromMult(multESD, fESDVtx[2] );
737  }
738  //
739  // Injection: it must come right after the normal reco since needs its results
740  if (GetDoInjection()) {
741  if (!fMultReco) InitMultReco(); // in principle, not needed, the reco is created above
742  fMultReco->SetRecType(AliITSMultRecBg::kBgInj);
743  fMultReco->Run(fRPTree, fESDVtx);
744  printf("Multiplicity from Injection:\n");
745  AliMultiplicity* mlt = fMultReco->GetMultiplicity();
746  if (mlt) mlt->Print();
747  hstat->Fill(kBinEntries + kEvProcInj + kEntriesPerBin*fCurrCentBin);
748  FillHistos(kBgInj,mlt);
749  }
750  //
751  // Rotation
752  if (GetDoRotation()) {
753  InitMultReco();
754  fMultReco->SetRecType(AliITSMultRecBg::kBgRot);
755  fMultReco->SetPhiRotationAngle(fPhiRot);
756  fMultReco->Run(fRPTree, fESDVtx);
757  printf("Multiplicity from Rotation:\n");
758  AliMultiplicity* mlt = fMultReco->GetMultiplicity();
759  if (mlt) mlt->Print();
760  hstat->Fill(kBinEntries + kEvProcRot + kEntriesPerBin*fCurrCentBin);
761  FillHistos(kBgRot,mlt);
762  }
763  //
764  // =============================================================================<<<
765  //
766  if (fMultReco) delete fMultReco;
767  fMultReco = 0;
768  //
769 }
770 
771 
772 //________________________________________________________________________
774 {
775  TH1* hstat;
776  TList *lst = dynamic_cast<TList*>(GetOutputData(1));
777  if (lst && (hstat=(TH1*)lst->FindObject("hStat"))) {
778  Info("Terminate","Registering used settings");
779  // fill used settings
780  hstat->Fill(kOneUnit,1.);
781  hstat->Fill(kCentVar,1);
782  hstat->Fill(kDPhi,fDPhiWindow);
783  hstat->Fill(kDTht,fDThetaWindow);
784  hstat->Fill(kNStd,fNStdDev);
785  hstat->Fill(kPhiShift,fDPhiShift);
786  hstat->Fill(kThtS2,fScaleDTBySin2T);
787  hstat->Fill(kThtCW,fCutOnDThetaX);
788  hstat->Fill(kPhiOvl,fPhiOverlapCut);
789  hstat->Fill(kZEtaOvl,fZetaOverlap);
790  hstat->Fill(kNoOvl,fRemoveOverlaps);
791  //
792  hstat->Fill(kPhiRot,fPhiRot);
793  hstat->Fill(kInjScl,fInjScale);
794  hstat->Fill(kEtaMin,fEtaMin);
795  hstat->Fill(kEtaMax,fEtaMax);
796  hstat->Fill(kZVMin,fZVertexMin);
797  hstat->Fill(kZVMax,fZVertexMax);
798  //
799  hstat->Fill(kDPiSCut,fDPhiSCut);
800  hstat->Fill(kNStdCut,fNStdCut);
801  hstat->Fill(kMCV0Scale, fMCV0Scale);
802  //
803  }
804  //
805 }
806 
807 
808 //________________________________________________________________________
810 {
811  Printf("Terminating...");
812  RegisterStat();
813 }
814 
815 
816 //_________________________________________________________________________
818 {
819  // create mult reconstructor
820  if (fMultReco) delete fMultReco;
821  fMultReco = new AliITSMultRecBg();
822  fMultReco->SetCreateClustersCopy(kTRUE);
823  fMultReco->SetScaleDThetaBySin2T(fScaleDTBySin2T);
824  fMultReco->SetNStdDev(fNStdDev);
825  fMultReco->SetPhiWindow( fDPhiWindow );
826  fMultReco->SetThetaWindow( fDThetaWindow );
827  fMultReco->SetPhiShift( fDPhiShift );
828  fMultReco->SetRemoveClustersFromOverlaps(fRemoveOverlaps);
829  fMultReco->SetPhiOverlapCut(fPhiOverlapCut);
830  fMultReco->SetZetaOverlapCut(fZetaOverlap);
831  fMultReco->SetHistOn(kFALSE);
832  fMultReco->SetRecType( AliITSMultRecBg::kData );
833 }
834 
835 //_________________________________________________________________________
837 {
838  // book custom histos, not related to specific tracklet type
839  TObjArray* histos = new TObjArray();
840  TH1F* hstat;
841  //
842  // ------------ job parameters, statistics ------------------------------>>>
844  hstat = new TH1F("hStat","Run statistics",nbs,0.5,nbs+0.5);
845  //
846  hstat->GetXaxis()->SetBinLabel(kEvTot0, "Ev.Tot0");
847  hstat->GetXaxis()->SetBinLabel(kEvTot, "Ev.Tot");
848  hstat->GetXaxis()->SetBinLabel(kOneUnit,"ScaleMerge");
849  hstat->GetXaxis()->SetBinLabel(kNWorkers,"Workers");
850  //
851  hstat->GetXaxis()->SetBinLabel(kCentVar, fUseCentralityVar.Data());
852  hstat->GetXaxis()->SetBinLabel(kDPhi, "#Delta#varphi");
853  hstat->GetXaxis()->SetBinLabel(kDTht, "#Delta#theta");
854  hstat->GetXaxis()->SetBinLabel(kNStd, "N.std");
855  hstat->GetXaxis()->SetBinLabel(kPhiShift,"#delta#varphi");
856  hstat->GetXaxis()->SetBinLabel(kThtS2,"scale #Delta#theta");
857  hstat->GetXaxis()->SetBinLabel(kPhiOvl,"#varpho_{Ovl}");
858  hstat->GetXaxis()->SetBinLabel(kZEtaOvl,"#z_{Ovl}");
859  hstat->GetXaxis()->SetBinLabel(kNoOvl, "rem.ovl");
860  //
861  hstat->GetXaxis()->SetBinLabel(kPhiRot,"#varphi_{rot}");
862  hstat->GetXaxis()->SetBinLabel(kInjScl,"inj");
863  hstat->GetXaxis()->SetBinLabel(kEtaMin,"#eta_{min}");
864  hstat->GetXaxis()->SetBinLabel(kEtaMax,"#eta_{max}");
865  hstat->GetXaxis()->SetBinLabel(kZVMin,"ZV_{min} cut");
866  hstat->GetXaxis()->SetBinLabel(kZVMax,"ZV_{max} cut");
867  //
868  hstat->GetXaxis()->SetBinLabel(kDPiSCut,"#Delta#varphi-#delta_{#phi} cut");
869  hstat->GetXaxis()->SetBinLabel(kNStdCut,"#Delta cut");
870  //
871  hstat->GetXaxis()->SetBinLabel(kMCV0Scale,"MC V0 scale");
872  //
873  hstat->GetXaxis()->SetBinLabel(kNEvSDMC,"MC PureSD");
874  //
875  for (int i=0;i<fNCentBins;i++) {
876  TString bnt = "b"; bnt+= i;
877  int offs = kBinEntries + i*kEntriesPerBin;
878  hstat->GetXaxis()->SetBinLabel(offs + kEvInMltBin, bnt+" Ev.In.CntBn");
879  hstat->GetXaxis()->SetBinLabel(offs + kEvProcData, bnt+" Ev.ProcData");
880  hstat->GetXaxis()->SetBinLabel(offs + kEvProcInj, bnt+" Ev.ProcInj");
881  hstat->GetXaxis()->SetBinLabel(offs + kEvProcRot, bnt+" Ev.ProcRot");
882  hstat->GetXaxis()->SetBinLabel(offs + kEvProcMix, bnt+" Ev.ProcMix");
883  hstat->GetXaxis()->SetBinLabel(offs + kEvCentBin, bnt+" Ev.CentBinNS");
884  hstat->GetXaxis()->SetBinLabel(offs + kEvPassPS, bnt+" Ev.CentBinPSTR");
885  hstat->GetXaxis()->SetBinLabel(offs + kEvPassVtx, bnt+" Ev.CentBinVtx");
886  //
887  }
888  //
889  hstat->Fill(kNWorkers);
890  //
891  AddHisto(histos,hstat,kHStat);
892  //
893  // ------------------------ events per centrality bin ----------------------
894  TH1F* hCentAx = new TH1F("EvCentr","Events per centrality",fNCentBins,fCentPerc.GetArray());
895  hCentAx->GetXaxis()->SetTitle("Centrality parameter");
896  AddHisto(histos,hCentAx,kHStatCent);
897  //
898  TH1F* hCentBin = new TH1F("EvCentrBin","Events per centrality bin",fNCentBins,-0.5,fNCentBins-0.5);
899  hCentBin->GetXaxis()->SetTitle("Centrality Bin");
900  AddHisto(histos,hCentBin,kHStatCentBin);
901  //
902  TH1F* hCentDistNoSel = new TH1F("EvCentrDistNoSel","Events per centrality Before selection",100,0,100);
903  hCentDistNoSel->GetXaxis()->SetTitle("Centrality percentile");
904  AddHisto(histos,hCentDistNoSel,kHCentDistNoSel);
905  //
906  TH1F* hCentDist = new TH1F("EvCentrDist","Events per centrality ",100,0,100);
907  hCentDist->GetXaxis()->SetTitle("Centrality percentile");
908  AddHisto(histos,hCentDist,kHCentDist);
909  //
910 
911  // ------------ job parameters, statistics ------------------------------<<<
912  //
913  // double etaMn=-3,etaMx=3;
914  double zMn=-30, zMx=30;
915  // int nEtaBins = TMath::Nint((etaMx-etaMn)/fEtaBinWidth);
916  // if (nEtaBins<1) nEtaBins = 1;
917  //
918  int nZVBins = int(zMx-zMn);
919  if (nZVBins<1) nZVBins = 1;
920  //
921  // Z vertex distribution for events before selection
922  TH1F* hzvns = new TH1F("zvNoSel","Z vertex before selection",nZVBins,zMn,zMx);
923  hzvns->GetXaxis()->SetTitle("Zvertex");
924  AddHisto(histos,hzvns,kHZVtxNoSel);
925  //
926  // V0 for events before selection
927  int nbmltV0 = 2500;
928  double maxmltV0 = 50000; //R+
929  //
930  TH1F* hnV0ns = new TH1F("V0NoSel","V0 signal Before Cent. Selection",nbmltV0,0,maxmltV0);
931  hnV0ns->GetXaxis()->SetTitle("V0 signal");
932  AddHisto(histos,hnV0ns,kHV0NoSel);
933  //
934  // N SPD2 clusters
935  int nbmltSPD2 = 2500;
936  double maxmltSPD2 = 25000; // R+
937  TH1F* hncl2ns = new TH1F("NClustersSPD2NoSel","N Clusters on SPD2 Before Cent Selection",nbmltSPD2,0,maxmltSPD2);
938  hncl2ns->GetXaxis()->SetTitle("N Clus SPD2");
939  AddHisto(histos,hncl2ns,kHNClSPD2NoSel);
940  //
941  int nbzdc=50;
942  double maxZDC=6000., maxZEM=2500.;
943  TH2F* hzdczemns = new TH2F("ZDCZEMNoSel","ZDC vs ZEM Before Cent Selection",
944  nbzdc,0,maxZEM,nbzdc,0,maxZDC);
945  hzdczemns->GetXaxis()->SetTitle("ZEM");
946  hzdczemns->GetXaxis()->SetTitle("ZDC");
947  AddHisto(histos,hzdczemns,kHZDCZEMNoSel);
948  //
949  // multiplicity distribution histos
950  TH1F* hmlt;
951  const int kMaxMlt = 10000;
952  hmlt = new TH1F("mltTrTPCITS","mltTPCITS",kMaxMlt,0,kMaxMlt);
953  AddHisto(histos,hmlt,kHMltEstTrITSTPC);
954  //
955  hmlt = new TH1F("mltTrITSSA","mltITSSA",kMaxMlt,0,kMaxMlt);
956  AddHisto(histos,hmlt,kHMltEstTrITSSA);
957  //
958  hmlt = new TH1F("mltTracklets","mltTracklets",kMaxMlt,0,kMaxMlt);
959  AddHisto(histos,hmlt,kHMltEstTr);
960  //
961  //
962  if (fUseMC) {
963  hmlt = new TH1F("mltMCTruth","mltMCTruth",kMaxMlt,0,kMaxMlt);
964  AddHisto(histos,hmlt,kHMltMC);
965  //
966  TH2F* hzvMCNoPS = new TH2F("zvMCNoPS","Z MC vertex Before Ph.Sel per Cent.Bin",nZVBins,zMn,zMx, fNCentBins, -0.5,fNCentBins-0.5);
967  hzvMCNoPS->GetXaxis()->SetTitle("ZvertexMC");
968  hzvMCNoPS->GetYaxis()->SetTitle("Cent.Bin ID");
969  AddHisto(histos,hzvMCNoPS,kHZVtxMCNoPhSel);
970  //
971  TH2F* hzvMCNoVS = new TH2F("zvMCNoVS","Z MC vertex Before Rec.Vtx.Selection per Cent.Bin",nZVBins,zMn,zMx, fNCentBins, -0.5,fNCentBins-0.5);
972  hzvMCNoVS->GetXaxis()->SetTitle("ZvertexMC");
973  hzvMCNoVS->GetYaxis()->SetTitle("Cent.Bin ID");
974  AddHisto(histos,hzvMCNoVS,kHZVtxMCNoVtSel);
975  //
976  TH2F* hzvMC = new TH2F("zvMC","Z vertex MC after Selection per Cent.Bin",nZVBins,zMn,zMx, fNCentBins, -0.5,fNCentBins-0.5);
977  hzvMC->GetXaxis()->SetTitle("ZvertexMC");
978  hzvMC->GetYaxis()->SetTitle("Cent.Bin ID");
979  AddHisto(histos,hzvMC,kHZVtxMC);
980  //
981  TH1* hNPrimEta2All = new TH1F("nPrimEta2All","NPrim |#eta|<2, allMC",1000,0.,1000.);
982  hNPrimEta2All->SetXTitle("NPrim");
983  AddHisto(histos,hNPrimEta2All,kHNPrimEta2All);
984  //
985  TH1* hNPrimEta2Vt = new TH1F("nPrimEta2Vt","NPrim |#eta|<2, VtxMC",1000,0.,1000.);
986  hNPrimEta2Vt->SetXTitle("NPrim");
987  AddHisto(histos,hNPrimEta2Vt,kHNPrimEta2Vt);
988  //
989  TH1* hNPrimEta2Sel = new TH1F("nPrimEta2Sel","NPrim |#eta|<2, Sel",1000,0.,1000.);
990  hNPrimEta2Sel->SetXTitle("NPrim");
991  AddHisto(histos,hNPrimEta2Sel,kHNPrimEta2Sel);
992  //
993  TH1* hNPrimEta2SelVt = new TH1F("nPrimEta2SelVt","NPrim |#eta|<2, Sel+VtxMC",1000,0.,1000.);
994  hNPrimEta2SelVt->SetXTitle("NPrim");
995  AddHisto(histos,hNPrimEta2SelVt,kHNPrimEta2SelVt);
996  //
997  }
998  //
999  TH2F* hcorrEta2 = new TH2F("nTrNprEta2","Ntr vs NPrim |#eta|<2, Sel+Vtx",100,0.,100.,100,0.,100.);
1000  hcorrEta2->SetXTitle("NPrim");
1001  hcorrEta2->SetYTitle("NTracklets");
1002  AddHisto(histos,hcorrEta2,kHNCorrMCEta2);
1003 
1004  TH2F* hzv = new TH2F("zv","Z vertex after Selection per Cent.Bin",nZVBins,zMn,zMx, fNCentBins, -0.5,fNCentBins-0.5);
1005  hzv->GetXaxis()->SetTitle("Zvertex");
1006  hzv->GetYaxis()->SetTitle("Cent.Bin ID");
1007  AddHisto(histos,hzv,kHZVtx);
1008  //
1009  // V0
1010  TH2F* hnV0 = new TH2F("V0","V0 signal per Cent.Bin ",nbmltV0,0,maxmltV0, fNCentBins, -0.5,fNCentBins-0.5);
1011  hnV0->GetXaxis()->SetTitle("V0 signal");
1012  hnV0->GetYaxis()->SetTitle("Cent.Bin ID");
1013  AddHisto(histos,hnV0,kHV0);
1014  //
1015  // N SPD2 clusters
1016  TH2F* hncl2 = new TH2F("NClustersSPD2","N Clusters on SPD2 per Cent.Bin",nbmltSPD2,0,maxmltSPD2, fNCentBins, -0.5,fNCentBins-0.5);
1017  hncl2->GetXaxis()->SetTitle("N Clus SPD2");
1018  hncl2->GetYaxis()->SetTitle("Cent.Bin ID");
1019  AddHisto(histos,hncl2,kHNClSPD2);
1020  //
1021  // ZDCZEM
1022  TH3F* hzdczem = new TH3F("ZDCZEM","ZDC vs ZEM per Cent.Bin",nbzdc,0,maxZEM,nbzdc,0,maxZDC, fNCentBins, -0.5,fNCentBins-0.5);
1023  hzdczem->GetXaxis()->SetTitle("ZEM");
1024  hzdczem->GetYaxis()->SetTitle("ZDC");
1025  AddHisto(histos,hzdczem,kHZDCZEM);
1026  //
1027  //----------------------------------------------------------------------
1028  TH2* hetaphi = new TH2F("etaphiTracklets","etaphiTracklets",50,-2.5,2.5, 200, 0,2*TMath::Pi());
1029  hetaphi->SetXTitle("#eta");
1030  hetaphi->SetYTitle("#phi");
1031  AddHisto(histos,hetaphi,kHEtaPhi);
1032  //----------------------------------------------------------------------
1033  int nEtaBinsS = TMath::Nint((fEtaMax-fEtaMin)/fEtaBinWidth);
1034  if (nEtaBinsS<1) nEtaBinsS = 1;
1035  //
1036  int nZVBinsS = int(fZVertexMax-fZVertexMin);
1037  if (nZVBinsS<1) nZVBinsS = 1;
1038 
1039  if (fUseMC) {
1040  // Z vertex vs Eta distribution for primaries
1041  char buffn[100],bufft[500];
1042  for (int ib=0;ib<fNCentBins;ib++) {
1043  sprintf(buffn,"b%d_zvEtaPrimMC",ib);
1044  sprintf(bufft,"bin%d Zvertex vs #eta PrimMC",ib);
1045  TH2F* hzvetap = new TH2F(buffn,bufft, nEtaBinsS,fEtaMin,fEtaMax,nZVBinsS,fZVertexMin,fZVertexMax);
1046  hzvetap->GetXaxis()->SetTitle("#eta");
1047  hzvetap->GetYaxis()->SetTitle("Zvertex");
1048  AddHisto(histos,hzvetap,kHZVEtaPrimMC+ib);
1049  //
1050  sprintf(buffn,"b%d_zvrecEtaPrimMCSel",ib);
1051  sprintf(bufft,"bin%d ZvertexR vs #eta PrimMC sel evs",ib);
1052  TH2F* hzvretap = new TH2F(buffn,bufft, nEtaBinsS,fEtaMin,fEtaMax,nZVBinsS,fZVertexMin,fZVertexMax);
1053  hzvretap->GetXaxis()->SetTitle("#eta");
1054  hzvretap->GetYaxis()->SetTitle("Zvertex");
1055  AddHisto(histos,hzvretap,kHZVrEtaPrimMC+ib);
1056  //
1057  sprintf(buffn,"b%d_dZV_ZVGen",ib);
1058  sprintf(bufft,"bin%d ZvRec-ZvGen vs ZvGen",ib);
1059  TH2F* hdzv = new TH2F(buffn,bufft,nZVBinsS,fZVertexMin,fZVertexMax, 200,-1.,1);
1060  hdzv->GetXaxis()->SetTitle("VZGen");
1061  hdzv->GetYaxis()->SetTitle("VZRec-VZGen");
1062  AddHisto(histos,hdzv,kHZVResMC+ib);
1063  //
1064  sprintf(buffn,"b%d_ptPrimMC",ib);
1065  sprintf(bufft,"bin%d pT PrimMC",ib);
1066  TH1F* hpt = new TH1F(buffn,bufft,1000,0.,5.);
1067  hpt->GetXaxis()->SetTitle("pT");
1068  hpt->GetYaxis()->SetTitle("counts");
1069  AddHisto(histos,hpt,kHPtPrimMC+ib);
1070  }
1071  //
1072  // <n> primaries according to MC generator
1073  TH1F* hnprimM = new TH1F("nPrimMean","<N> primaries",fNCentBins, -0.5,fNCentBins-0.5);
1074  hnprimM->GetXaxis()->SetTitle("Cent.Bin ID");
1075  AddHisto(histos,hnprimM,kHNPrimMeanMC);
1076  //
1077  // <n> primaries per part.pair according to MC generator
1078  TH1F* hnprim2partM = new TH1F("nPrim2Part","<N> primaries per part.pair",fNCentBins, -0.5,fNCentBins-0.5);
1079  hnprim2partM->GetXaxis()->SetTitle("Cent.Bin ID");
1080  AddHisto(histos,hnprim2partM,kHNPrim2PartMC);
1081  //
1082  // <n> primaries per part.pair vs npart.pair according to MC generator
1083  TH2F* hnprim2partNp = new TH2F("nPrim2Part_vs_NPart","<N> primaries per part.pair vs N part.pairs",105,0,210,200,0,40);
1084  hnprim2partNp->GetXaxis()->SetTitle("N.part.pairs");
1085  hnprim2partNp->GetYaxis()->SetTitle("N.prim/N.part.pairs");
1086  AddHisto(histos,hnprim2partNp,kHNPrim2PartNpMC);
1087  //
1088  // <n> primaries per b.coll vs npart.pair according to MC generator
1089  TH2F* hnprim2BCollNp = new TH2F("nPrim2BColl_vs_NPart","<N> primaries per bin.coll vs N part.pairs",105,0,210,200,0,40);
1090  hnprim2BCollNp->GetXaxis()->SetTitle("N.part.pairs");
1091  hnprim2BCollNp->GetYaxis()->SetTitle("N.prim/N.bin.coll.");
1092  AddHisto(histos,hnprim2BCollNp,kHNPrim2BCollNpMC);
1093  //
1094  // <n> primaries per bin.coll. according to MC generator
1095  TH1F* hnprim2BCollM = new TH1F("nPrim2BColl","<N> primaries per bin.coll",fNCentBins, -0.5,fNCentBins-0.5);
1096  hnprim2BCollM->GetXaxis()->SetTitle("Cent.Bin ID");
1097  AddHisto(histos,hnprim2BCollM,kHNPrim2BCollMC);
1098  //
1099  // n participants according to MC generator
1100  TH2F* hnpart = new TH2F("nPart","N participant pairs",210,0,210,fNCentBins, -0.5,fNCentBins-0.5);
1101  hnpart->GetXaxis()->SetTitle("N part. pairs");
1102  hnpart->GetYaxis()->SetTitle("Cent.Bin ID");
1103  AddHisto(histos,hnpart,kHNPartMC);
1104  //
1105  // <n> participants according to MC generator
1106  TH1F* hnpartM = new TH1F("nPartMean","<N> participant pairs",fNCentBins, -0.5,fNCentBins-0.5);
1107  hnpartM->GetXaxis()->SetTitle("Cent.Bin ID");
1108  AddHisto(histos,hnpartM,kHNPartMeanMC);
1109  //
1110  // n bin coll. according to MC generator
1111  TH2F* hnbcoll = new TH2F("nBColl","N bin. coll",2000,0,2000,fNCentBins, -0.5,fNCentBins-0.5);
1112  hnbcoll->GetXaxis()->SetTitle("N bin. coll");
1113  hnbcoll->GetYaxis()->SetTitle("Cent.Bin ID");
1114  AddHisto(histos,hnbcoll,kHNBCollMC);
1115  //
1116  // <n> bin col according to MC generator
1117  TH1F* hnbcollM = new TH1F("nBCollMean","<N> bin.colls",fNCentBins, -0.5,fNCentBins-0.5);
1118  hnbcollM->GetXaxis()->SetTitle("Cent.Bin ID");
1119  AddHisto(histos,hnbcollM,kHNBCollMeanMC);
1120  //
1121  }
1122  //
1123  // --------------------------------------------------
1124  if (fUseMC) {
1125  int npdg = sizeof(fgkPDGNames)/sizeof(char*);
1126  TH2F* hpdgP = new TH2F("pdgPrim","primary PDG",npdg,0,npdg,fNCentBins, -0.5,fNCentBins-0.5);
1127  AddHisto(histos,hpdgP,kHPrimPDG);
1128  TH2F* hpdgS = new TH2F("pdgSec","secondary PDG",npdg,0,npdg,fNCentBins, -0.5,fNCentBins-0.5);
1129  AddHisto(histos,hpdgS,kHSecPDG);
1130  TH2F* hpdgPP = new TH2F("pdgPrimPar","primary parent PDG ",npdg,0,npdg,fNCentBins, -0.5,fNCentBins-0.5);
1131  AddHisto(histos,hpdgPP,kHPrimParPDG);
1132  TH2F* hpdgSP = new TH2F("pdgSecPar","secondary parent PDG",npdg,0,npdg,fNCentBins, -0.5,fNCentBins-0.5);
1133  AddHisto(histos,hpdgSP,kHSecParPDG);
1134  TH2F* hpdgMC = new TH2F("pdgPrimMC","primary PDG",npdg,0,npdg,fNCentBins, -0.5,fNCentBins-0.5);
1135  AddHisto(histos,hpdgMC,kHPdgPrimMC);
1136  for (int i=0;i<npdg;i++) {
1137  hpdgP->GetXaxis()->SetBinLabel(i+1,fgkPDGNames[i]);
1138  hpdgS->GetXaxis()->SetBinLabel(i+1,fgkPDGNames[i]);
1139  hpdgPP->GetXaxis()->SetBinLabel(i+1,fgkPDGNames[i]);
1140  hpdgSP->GetXaxis()->SetBinLabel(i+1,fgkPDGNames[i]);
1141  hpdgMC->GetXaxis()->SetBinLabel(i+1,fgkPDGNames[i]);
1142  }
1143  }
1144  //
1145  // -------------------------------------------------
1146  TH2F* hclinf=0;
1147  hclinf = new TH2F("cl0InfoUsed","#phi vs Z of used clusters, Lr0",64,-16,16, 80,0,2*TMath::Pi());
1148  AddHisto(histos,hclinf,kHClUsedInfoL0);
1149  hclinf = new TH2F("cl1InfoUsed","#phi vs Z of used clusters, Lr1",64,-16,16, 2*80,0,2*TMath::Pi());
1150  AddHisto(histos,hclinf,kHClUsedInfoL1);
1151  hclinf = new TH2F("cl0InfoAll","#phi vs Z of all clusters, Lr0",64,-16,16, 80,0,2*TMath::Pi());
1152  AddHisto(histos,hclinf,kHClAllInfoL0);
1153  hclinf = new TH2F("cl1InfoAll","#phi vs Z of all clusters, Lr1",64,-16,16, 2*80,0,2*TMath::Pi());
1154  AddHisto(histos,hclinf,kHClAllInfoL1);
1155  //
1156  // -------------------------------------------------
1157  histos->SetOwner(kFALSE);
1158  //
1159  return histos;
1160 }
1161 
1162 //_________________________________________________________________________
1164 {
1165  // book standard set of histos attaching the pref in front of the name/title
1166  //
1167  const int kNDPhiBins = 100;
1168  const int kNDThtBins = 100;
1169  int nDistBins = int(fNStdDev)*5;
1170  //
1171  int nEtaBins = TMath::Nint((fEtaMax-fEtaMin)/fEtaBinWidth);
1172  if (nEtaBins<1) nEtaBins = 1;
1173  //
1174  int nZVBins = int(fZVertexMax-fZVertexMin);
1175  if (nZVBins<1) nZVBins = 1;
1176  float dphir = fDPhiWindow*TMath::Sqrt(fNStdDev);
1177  float dthtr = fDThetaWindow*TMath::Sqrt(fNStdDev);
1178  //
1179  TObjArray* histos = new TObjArray();
1180  TH2F* h2;
1181  TH1F* h1;
1182  char buffn[100],bufft[500];
1183  //
1184  for (int ib=0;ib<fNCentBins;ib++) {
1185  //
1186  int offs = ib*kNStandardH;
1187  if (selHistos & (0x1<<kHEtaZvCut) ) {
1188  sprintf(buffn,"b%d_%s_ZvEtaCutT",ib,pref);
1189  sprintf(bufft,"bin%d (%s) Zv vs Eta with tracklet cut",ib,pref);
1190  h2 = new TH2F(buffn,bufft,nEtaBins,fEtaMin,fEtaMax, nZVBins, fZVertexMin,fZVertexMax);
1191  h2->GetXaxis()->SetTitle("#eta");
1192  h2->GetYaxis()->SetTitle("Zv");
1193  AddHisto(histos,h2,offs+kHEtaZvCut);
1194  }
1195  //
1196  if (selHistos & (0x1<<kHDPhiDTheta) ) {
1197  sprintf(buffn,"b%d_%s_dPhidTheta",ib,pref);
1198  sprintf(bufft,"bin%d (%s) #Delta#theta vs #Delta#varphi",ib,pref);
1199  h2 = new TH2F(buffn,bufft,kNDPhiBins,-dphir,dphir,kNDThtBins,-dthtr,dthtr);
1200  h2->GetXaxis()->SetTitle("#Delta#varphi [rad]");
1201  h2->GetYaxis()->SetTitle("#Delta#theta [rad]");
1202  AddHisto(histos,h2,offs+kHDPhiDTheta);
1203  }
1204  //
1205  if (selHistos & (0x1<<kHDPhiSDThetaX) ) {
1206  sprintf(buffn,"b%d_%s_dPhiSdThetaX",ib,pref);
1207  sprintf(bufft,"bin%d (%s) #Delta#theta%s vs #Delta#varphi-#delta_{#varphi}",ib,pref,fScaleDTBySin2T ? "/sin^{2}(#theta)":"");
1208  h2 = new TH2F(buffn,bufft,kNDPhiBins,-dphir,dphir,kNDThtBins,-dthtr,dthtr);
1209  h2->GetXaxis()->SetTitle("#Delta#varphi-#delta_{#varphi} [rad]");
1210  sprintf(bufft,"#Delta#theta%s",fScaleDTBySin2T ? "/sin^{2}(#theta)":"");
1211  h2->GetYaxis()->SetTitle(bufft);
1212  AddHisto(histos,h2,offs+kHDPhiSDThetaX);
1213  }
1214  //
1215  if (selHistos & (0x1<<kHWDist) ) {
1216  sprintf(buffn,"b%d_%s_WDist",ib,pref);
1217  sprintf(bufft,"bin%d #Delta=[(#Delta#varphi-#delta_{#varphi})/#sigma#varphi]^{2}+"
1218  "[#Delta#theta%s/#sigma#theta]^{2}",ib,fScaleDTBySin2T ? "*sin^{-2}(#theta)":"");
1219  h1 = new TH1F(buffn,bufft,nDistBins,0,fNStdDev);
1220  sprintf(bufft,"#Delta=[(#Delta#varphi-#delta_{#varphi})/#sigma#varphi]^{2}+"
1221  "[#Delta#theta%s/#sigma#theta]^{2}",fScaleDTBySin2T ? "*sin^{-2}(#theta)":"");
1222  h1->GetXaxis()->SetTitle(bufft);
1223  AddHisto(histos,h1,offs+kHWDist);
1224  }
1225  //
1226  /*
1227  if (selHistos & (0x1<<kHEtaZvSPD1) ) {
1228  sprintf(buffn,"b%d_%s_ZvEtaSPD1",ib,pref);
1229  sprintf(bufft,"bin%d (%s) Zv vs Eta SPD1 clusters",ib,pref);
1230  h2 = new TH2F(buffn,bufft,nEtaBins,fEtaMin,fEtaMax, nZVBins, fZVertexMin,fZVertexMax);
1231  h2->GetXaxis()->SetTitle("#eta");
1232  h2->GetYaxis()->SetTitle("Zv");
1233  AddHisto(histos,h2,offs+kHEtaZvSPD1);
1234  }
1235  */
1236  //
1237  if (selHistos & (0x1<<kHWDvEta) ) {
1238  sprintf(buffn,"b%d_%s_WDvsEta",ib,pref);
1239  sprintf(bufft,"bin%d (%s) Wdist vs Eta",ib,pref);
1240  h2 = new TH2F(buffn,bufft,nEtaBins,fEtaMin,fEtaMax, 2*nDistBins,0,fNStdDev);
1241  h2->GetXaxis()->SetTitle("#eta");
1242  h2->GetYaxis()->SetTitle("Wdist");
1243  AddHisto(histos,h2,offs+kHWDvEta);
1244  }
1245  //
1246  }
1247  //
1248  histos->SetOwner(kFALSE);
1249  return histos;
1250 }
1251 
1252 //_________________________________________________________________________
1254 {
1255  // add single histo to the set
1256  if (at>=0) histos->AddAtAndExpand(h,at);
1257  else histos->Add(h);
1258  fOutput->Add(h);
1259 }
1260 
1261 //_________________________________________________________________________
1262 void AliTrackletTaskMulti::FillHistos(Int_t type, const AliMultiplicity* mlt)
1263 {
1264  // fill histos of given type
1265  if (!mlt) return;
1266  //
1267  TObjArray* histos = 0;
1268  if (type == kData) histos = fHistosTrData;
1269  else if (type == kBgInj) histos = fHistosTrInj;
1270  else if (type == kBgRot) histos = fHistosTrRot;
1271  //
1272  Bool_t fillMC = (type==kData) && fUseMC && fStack;
1273  //
1274  //
1275  //---------------------------------------- CHECK ------------------------------>>>
1276  TArrayF vtxMC;
1277  AliGenHijingEventHeader* pyHeader = 0;
1278  //
1279  if (fUseMC) {
1280  pyHeader = (AliGenHijingEventHeader*) fMCEvent->GenEventHeader();//header->GenEventHeader();
1281  pyHeader->PrimaryVertex(vtxMC);
1282  }
1283  //---------------------------------------- CHECK ------------------------------<<<
1284  //
1285  if (!histos) return;
1286  int ntr = mlt->GetNumberOfTracklets();
1287  fNTreta2 = 0;
1288  for (int itr=ntr;itr--;) {
1289  //
1290  //--------------------------------------- REWEIGHT ---------------------------->>>
1291  // R+HACK
1292  //
1293  double weight = 1.;
1294  if (fUseMC) {
1295  if (mlt->GetLabel(itr,0) == mlt->GetLabel(itr,1)) {
1296  int lab = mlt->GetLabel(itr,0);
1297  if (lab < 0) continue;
1298  TParticle* part = fStack->Particle(lab);
1299  if (!part) continue;
1300  weight = part->GetWeight();
1301  }
1302  else {
1303  for (Int_t ilay = 0; ilay < 2; ilay++) {
1304  int lab = mlt->GetLabel(itr,ilay);
1305  if (lab < 0) continue;
1306  TParticle* part = fStack->Particle(lab);
1307  if (!part) continue;
1308  weight *= part->GetWeight();
1309  }
1310  }
1311  }
1312  //---------------------------------------- CHECK ------------------------------>>>
1313  /*
1314  if (fUseMC) {
1315  Bool_t reject = kFALSE;
1316  while(1) {
1317  int lab0 = mlt->GetLabel(itr,0);
1318  int lab1 = mlt->GetLabel(itr,1);
1319  if (lab0!=lab1) break;
1320  if (!fStack->IsPhysicalPrimary(lab0)) break;
1321  //
1322  TParticle* part = fStack->Particle(lab0);
1323  Float_t dz = part->Vz() - vtxMC[2];
1324  if (TMath::Abs(dz)<1e-6) break;
1325  reject = kTRUE;
1326  break;
1327  }
1328  if (reject) continue;
1329  }
1330  */
1331  //---------------------------------------- CHECK ------------------------------<<<
1332  //
1333  double theta = mlt->GetTheta(itr);
1334  double phi = mlt->GetPhi(itr); // R+HACK
1335  double eta = -TMath::Log(TMath::Tan(theta/2));
1336  //
1337  double dtheta = mlt->GetDeltaTheta(itr);
1338  double dThetaX = dtheta;
1339  if (fScaleDTBySin2T) {
1340  double sint = TMath::Sin(theta);
1341  dThetaX /= (sint*sint);
1342  }
1343  if (fCutOnDThetaX && TMath::Abs(dThetaX)>fDThetaWindow) continue;
1344  //
1345  // double phi = mlt->GetPhi(itr);
1346  double dphi = mlt->GetDeltaPhi(itr);
1347  double dist = mlt->CalcDist(itr);
1348  double dphiS = TMath::Abs(dphi) - fDPhiShift; if (dphi<0) dphiS = -dphiS;
1349  //
1350  if (dist<fNStdCut && dphiS<fDPhiSCut && TMath::Abs(eta)<2) fNTreta2++;
1351  //
1352  //
1353  if (eta<fEtaMin || eta>fEtaMax) continue;
1354  if (phi<fPhiMin || phi>fPhiMax) continue; // R+HACK
1355  //
1356  Bool_t isSig = FillHistosSet(histos,eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist,weight); // R+HACK
1357  //
1358  if (type==kData && isSig) {
1359  ((TH2*)fHistosCustom->UncheckedAt(kHEtaPhi))->Fill(eta,mlt->GetPhi(itr),weight); // R+HACK
1360  }
1361  //
1362  // special handling for mc info
1363  if (fillMC && fStack) {
1364  int lab0 = mlt->GetLabel(itr,0);
1365  int lab1 = mlt->GetLabel(itr,1);
1366  int typeMC = 2; // comb.bg.
1367  if (lab0 == lab1) typeMC = fStack->IsPhysicalPrimary(lab0) ? 0:1; // prim or sec
1368  if (typeMC==0) FillHistosSet(fHistosTrPrim,eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist,weight); // primary // R+HACK
1369  else if (typeMC==1) FillHistosSet(fHistosTrSec, eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist,weight); // secondary // R+HACK
1370  else {
1371  FillHistosSet(fHistosTrComb,eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist,weight); // comb // R+HACK
1372  // for combinatorals fill also the uncorrelated part
1373  if (fMultReco) {
1374  float *trl = fMultReco->GetTracklet(itr);
1375  int clId0 = (int)trl[AliITSMultReconstructor::kClID1];
1376  int clId1 = (int)trl[AliITSMultReconstructor::kClID2];
1377  float *clLabs0 = fMultReco->GetClusterOfLayer(0,clId0) + AliITSMultReconstructor::kClMC0;
1378  float *clLabs1 = fMultReco->GetClusterOfLayer(1,clId1) + AliITSMultReconstructor::kClMC0;
1379  if (!HaveCommonParent(clLabs0,clLabs1))
1380  FillHistosSet(fHistosTrCombU,eta,/*phi,theta,*/dphi,dtheta,dThetaX,dist,weight); // R+HACK
1381  }
1382  } // combinatorials
1383 
1384  if (dist<fNStdCut) {
1385  if (dphiS<fDPhiSCut) FillSpecies(typeMC,lab0);
1386  }
1388  }
1389  }
1390  //
1391  if (type==kData) {
1392  TH2* hcorr = ((TH2*)fHistosCustom->UncheckedAt(kHNCorrMCEta2));
1393  if (hcorr) hcorr->Fill(fNPrimMCeta2,fNTreta2);
1394  }
1395  //
1396  //-------------------------------------------------------------TMP RS - singles ------->>>
1397  /*
1398  int offsH = fCurrCentBin*kNStandardH;
1399  TH2* hSingles = (TH2*)histos->UncheckedAt(offsH+kHEtaZvSPD1);
1400  if (hSingles) {
1401  int nclS = mlt->GetNumberOfSingleClusters();
1402  double *thtS = mlt->GetThetaSingle();
1403  for (int ics=nclS;ics--;) {
1404  double etaS = -TMath::Log(TMath::Tan(thtS[ics]/2));
1405  if (etaS<fEtaMin || etaS>fEtaMax) continue;
1406  hSingles->Fill(etaS,fESDVtx[2]);
1407  }
1408  }
1409  */
1410  //-------------------------------------------------------------TMP RS - singles -------<<<
1411  //
1412 }
1413 
1414 //_________________________________________________________________________
1416 {
1417  // fill all MC primaries Zv vs Eta
1418  if (!fUseMC || !fStack || !fMCEvent) return;
1419  //
1420  float zv = fVtxMC[2]; //fVtxOK ? fESDVtx[2] : fVtxMC[2];
1421  float zvr = fESDVtx[2];
1422  /*
1423  const double kSafeZv = 1e-3
1424  if (zv<fZVertexMin+kSafeZv) zv = fZVertexMin+kSafeZv;
1425  if (zv>fZVertexMax-kSafeZv) zv = fZVertexMax-kSafeZv;
1426  */
1427  //
1428  int ntr = fStack->GetNtrack();
1429  TH2* hprimEtaZ = (TH2F*)fHistosCustom->UncheckedAt(kHZVEtaPrimMC+fCurrCentBin);
1430  TH2* hprimEtaZr = (TH2F*)fHistosCustom->UncheckedAt(kHZVrEtaPrimMC+fCurrCentBin);
1431  TH1* hprimPt = (TH1F*)fHistosCustom->UncheckedAt(kHPtPrimMC+fCurrCentBin); // R+HACK
1432  TH2 *hprimPdg = (TH2F*)fHistosCustom->UncheckedAt(kHPdgPrimMC); // R+HACK
1433  if (fVtxOK) {
1434  TH2* hdzv = (TH2F*)fHistosCustom->UncheckedAt(kHZVResMC+fCurrCentBin);
1435  hdzv->Fill(fVtxMC[2],fESDVtx[2]-fVtxMC[2]);
1436  }
1437  int nprim = 0;
1438  fNPrimMCeta2 = 0;
1439  float nprimSel = 0;
1440  for (int itr=ntr;itr--;) {
1441  if (!fStack->IsPhysicalPrimary(itr)) continue;
1442  AliMCParticle *mcpart = (AliMCParticle*)fMCEvent->GetTrack(itr);
1443  Float_t theta = mcpart->Theta();
1444  if (theta<1e-6 || theta>TMath::Pi()-1e-6) continue;
1445  Float_t eta = mcpart->Eta();
1446  Float_t pt = mcpart->Pt();
1447  //
1448  //--------------------------------------- REWEIGHT ---------------------------->>>
1449  // R+HACK
1450  //
1451  TParticle* part = fStack->Particle(itr);
1452  if (!part) continue;
1453  double weight = part->GetWeight();
1454  int pdgCode = TMath::Abs(part->GetPdgCode());
1455  int pdgBin = GetPdgBin(pdgCode);
1456  //
1457  //--------------------------------------- MC INFO ---------------------------->>>
1458  // R+HACK
1459  //
1460  // inclusive primary charged spectrum
1461  if (mcpart->Charge() && TMath::Abs(eta) < 1.)
1462  hprimPt->Fill(pt,weight);
1463  // inclusive primary particle composition
1464  if (TMath::Abs(eta) < 1.)
1465  hprimPdg->Fill(pdgBin,fCurrCentBin,weight); // R+HACK
1466  //
1467  if (!mcpart->Charge()) continue;
1468  //
1469  //---------------------------------------- CHECK ------------------------------>>>
1470  /*
1471  Float_t dz = part->Zv() - vtxMC[2];
1472  if (TMath::Abs(dz)>1e-6) continue; // reject
1473  */
1474  //---------------------------------------- CHECK ------------------------------<<<
1475  //
1476  // if (eta<fEtaMin || eta>fEtaMax) continue;
1477  if (TMath::Abs(eta)<2) fNPrimMCeta2++;
1478  hprimEtaZ->Fill(eta, zv, weight);
1479  if (fVtxOK && fIsSelected) {
1480  hprimEtaZr->Fill(eta, zvr, weight);
1481  if (eta>fEtaMin && eta<fEtaMax) nprimSel++;
1482  }
1483  nprim++;
1484  }
1485  ((TH1F*)fHistosCustom->UncheckedAt(kHNPrimEta2All))->Fill(fNPrimMCeta2);
1486  if (fVtxOK) ((TH1F*)fHistosCustom->UncheckedAt(kHNPrimEta2Vt))->Fill(fNPrimMCeta2);
1487  if (fIsSelected) ((TH1F*)fHistosCustom->UncheckedAt(kHNPrimEta2Sel))->Fill(fNPrimMCeta2);
1488  if (fVtxOK && fIsSelected) ((TH1F*)fHistosCustom->UncheckedAt(kHNPrimEta2SelVt))->Fill(fNPrimMCeta2);
1489  //
1490  if (fIsSelected) ((TH1F*)fHistosCustom->UncheckedAt(kHMltMC))->Fill(nprimSel);
1491  //
1492  ((TH1*)fHistosCustom->UncheckedAt(kHNPrimMeanMC))->Fill(fCurrCentBin,nprim);
1493  if (fNPart>0) {
1494  ((TH1*)fHistosCustom->UncheckedAt(kHNPrim2PartMC))->Fill(fCurrCentBin,nprim/fNPart);
1495  ((TH2*)fHistosCustom->UncheckedAt(kHNPrim2PartNpMC))->Fill(fNPart,nprim/fNPart);
1496  ((TH2*)fHistosCustom->UncheckedAt(kHNPartMC))->Fill(fNPart,fCurrCentBin);
1497  ((TH1*)fHistosCustom->UncheckedAt(kHNPartMeanMC))->Fill(fCurrCentBin,fNPart);
1498  }
1499  if (fNBColl>0) {
1500  ((TH1*)fHistosCustom->UncheckedAt(kHNPrim2BCollMC))->Fill(fCurrCentBin,nprim/fNBColl);
1501  ((TH2*)fHistosCustom->UncheckedAt(kHNPrim2BCollNpMC))->Fill(fNPart,nprim/fNBColl);
1502  ((TH2*)fHistosCustom->UncheckedAt(kHNBCollMC))->Fill(fNBColl,fCurrCentBin);
1503  ((TH1*)fHistosCustom->UncheckedAt(kHNBCollMeanMC))->Fill(fCurrCentBin,fNBColl);
1504  }
1505  //
1506 }
1507 
1508 //_________________________________________________________________________
1510  //double /*phi*/,double /*theta*/,
1511  double dphi,double dtheta,double dThetaX,
1512  double dist,double weight)
1513 {
1514  // fill standard set of histos
1515  Bool_t res = kFALSE;
1516  //
1517  int offs = fCurrCentBin*kNStandardH;
1518  //
1519  if (histos->UncheckedAt(kHWDvEta))
1520  ((TH2*)histos->UncheckedAt(offs+kHWDvEta))->Fill(eta,dist,weight);
1521  //
1522  if (dist>fNStdDev) return res;
1523  //
1524  double dphiS = TMath::Abs(dphi) - fDPhiShift;
1525  if (dphi<0) dphiS = -dphiS;
1526  //
1527  if (histos->UncheckedAt(offs+kHDPhiSDThetaX))
1528  ((TH2*)histos->UncheckedAt(offs+kHDPhiSDThetaX))->Fill(dphiS,dThetaX,weight);
1529  //
1530  if (histos->UncheckedAt(offs+kHDPhiDTheta))
1531  ((TH2*)histos->UncheckedAt(offs+kHDPhiDTheta))->Fill(dphi,dtheta,weight);
1532  //
1533  if (histos->UncheckedAt(kHWDist))
1534  ((TH2*)histos->UncheckedAt(offs+kHWDist))->Fill(dist,weight);
1535  //
1536  if (dist<fNStdCut && dphiS<fDPhiSCut && histos->UncheckedAt(offs+kHEtaZvCut)) {
1537  ((TH2*)histos->UncheckedAt(offs+kHEtaZvCut))->Fill(eta,fESDVtx[2],weight);
1538  res = kTRUE;
1539  }
1540  //
1541  return res;
1542 }
1543 
1544 //__________________________________________________________________
1546 {
1547  // fill PDGcode
1548  TH2 *hPart=0,*hParent=0; // R+HACK (it was TH1 *)
1549  if (primsec==0) {
1550  hPart = (TH2*)fHistosCustom->UncheckedAt(kHPrimPDG); // R+HACK (it was TH1 *)
1551  hParent = (TH2*)fHistosCustom->UncheckedAt(kHPrimParPDG); // R+HACK (it was TH1 *)
1552  }
1553  else if (primsec==1) {
1554  hPart = (TH2*)fHistosCustom->UncheckedAt(kHSecPDG); // R+HACK (it was TH1 *)
1555  hParent = (TH2*)fHistosCustom->UncheckedAt(kHSecParPDG); // R+HACK (it was TH1 *)
1556  }
1557  else return;
1558  int ntr = fStack->GetNtrack();
1559  TParticle* part = fStack->Particle(id);
1560  double weight = part->GetWeight();
1561  int pdgCode = TMath::Abs(part->GetPdgCode());
1562  int pdgBin = GetPdgBin(pdgCode);
1563  int parID = part->GetFirstMother();
1564  int pdgCodePar = -1;
1565  int pdgBinPar = -1;
1566  while (parID>=0 && parID<ntr) {
1567  part = fStack->Particle(parID);
1568  pdgCodePar = TMath::Abs(part->GetPdgCode());
1569  parID = part->GetFirstMother();
1570  }
1571  if (pdgCodePar>0) pdgBinPar = GetPdgBin(pdgCodePar);
1572  //
1573  hPart->Fill(pdgBin,fCurrCentBin,weight); // R+HACK
1574  hParent->Fill(pdgBinPar,fCurrCentBin,weight); // R+HACK
1575  //
1576 }
1577 
1578 //_________________________________________________________________________
1580 {
1581  // calculate centrality bin
1582  //R+HACK for (int i=0;i<fNCentBins;i++) if (perc>=fCentPerc[i] && perc<=fCentPerc[i+1]) return i;
1583  for (int i=0;i<fNCentBins;i++) if (perc>=fCentPerc[i] && perc<fCentPerc[i+1]) return i; // R+FIX
1584  return -1;
1585 }
1586 
1587 //_________________________________________________________________________
1589 {
1590  // return my pdg bin
1591  int ncodes = sizeof(fgkPDGCodes)/sizeof(int);
1592  int pdgBin=0;
1593  for (pdgBin=0;pdgBin<ncodes;pdgBin++) if (pdgCode==fgkPDGCodes[pdgBin]) break;
1594  if (pdgBin>=ncodes) {
1595  if (float(pdgCode)>1e9) pdgBin = ncodes; // nuclei
1596  else pdgBin = ncodes+1; // unknown
1597  }
1598  return pdgBin;
1599 }
1600 
1601 //_________________________________________________________________________
1602 Bool_t AliTrackletTaskMulti::HaveCommonParent(const float* clLabs0,const float* clLabs1)
1603 {
1604  // do 2 clusters have common parrent
1605  const int kMaxPar = 50;
1606  static int pars[2][50];
1607  int npars[2]={0,0};
1608  const float *labs[2] = {clLabs0,clLabs1};
1609  int ntr = fStack->GetNtrack();
1610  for (int il=0;il<2;il++) {
1611  for (int ilb=0;ilb<3;ilb++) {
1612  int lbl = (int)labs[il][ilb];
1613  if (lbl<0 || lbl>=ntr) continue;
1614  //
1615  while (npars[il]<kMaxPar-1) {
1616  pars[il][ npars[il]++ ] = lbl;
1617  TParticle* part = fStack->Particle(lbl);
1618  if (!part) break;
1619  lbl = part->GetFirstMother();
1620  if (lbl<1 || lbl>=ntr) break;
1621  }
1622  }
1623  }
1624  // compare array of parents
1625  for (int i0=npars[0];i0--;) for (int i1=npars[1];i1--;) if (pars[0][i0]==pars[1][i1]) return kTRUE;
1626  return kFALSE;
1627 }
1628 
1629 
1630 //_________________________________________________________________________
1632 {
1633  // R+DONTFORGET
1634  // fill reconstructable tracklets hitsos
1635  static TArrayI trInd;
1636  static TArrayF trWeight;
1637  static TBits isPrimArr;
1638  //
1639  if (!fMultReco || !fMultReco->IsRecoDone()) {AliInfo("To check reconstructables the reco had to be requested"); return;}
1640  if (!fStack) {AliInfo("MC Stack is not availalble"); return;}
1641  const double kPtMin = 0.05;
1642  //
1643  TClonesArray *clArr[2];
1644  for (int ilr=0;ilr<2;ilr++) {
1645  clArr[ilr] = fMultReco->GetClustersOfLayer(ilr);
1646  if (!clArr[ilr]) {AliInfo("Clusters are not available"); return;}
1647  }
1648  //
1649  int ntr = fStack->GetNtrack();
1650  if (!ntr) return;
1651  trInd.Reset();
1652  if (trInd.GetSize()<ntr) trInd.Set(ntr);
1653  isPrimArr.ResetAllBits();
1654  // count track wich may be reconstructable
1655  //
1656  int ntrStore=0,ntrStorePrim=0;
1657  Int_t *trIndArr = trInd.GetArray();
1658  Float_t *trWeightArr = trWeight.GetArray();
1659  for (Int_t it=0; it<ntr; it++) {
1660  TParticle* part = fStack->Particle(it);
1661  if (TMath::Abs(part->Eta())>2.2) continue;
1662  if (TMath::Abs(part->Pt())<kPtMin) continue;
1663  if (fStack->IsPhysicalPrimary(it)) {
1664  isPrimArr.SetBitNumber(it);
1665  ntrStorePrim++;
1666  }
1667  else { // check if secondary is worth cheking
1668  TParticlePDG* pdgPart = part->GetPDG();
1669  if (TMath::Abs(pdgPart->Charge())!=3) continue;
1670  if (part->R()>5.) continue;
1671  }
1672  trIndArr[it] = ++ntrStore;
1673  trWeightArr[it] = part->GetWeight();
1674  }
1675  //
1676  AliInfo(Form("Selected %d MC particles (%d prim) out of %d in the stack\n",ntrStore,ntrStorePrim,ntr));
1677  //
1678  const int kMaxCl=3;
1679  AliITSRecPoint **clIndL[2];
1680  clIndL[0] = new AliITSRecPoint*[kMaxCl*ntrStore]; // max 2 clusters per layer
1681  clIndL[1] = new AliITSRecPoint*[kMaxCl*ntrStore]; // max 2 clusters per layer
1682  memset(clIndL[0],0,kMaxCl*ntrStore*sizeof(AliITSRecPoint*));
1683  memset(clIndL[1],0,kMaxCl*ntrStore*sizeof(AliITSRecPoint*));
1684  //
1685  for (int ilr=0;ilr<2;ilr++) {
1686  TClonesArray *clusters = clArr[ilr];
1687  int ncl = clusters->GetEntriesFast();
1688  for (int icl=ncl;icl--;) {
1689  AliITSRecPoint *cl = (AliITSRecPoint*)clusters->UncheckedAt(icl);
1690  for (int ilb=3;ilb--;) {
1691  int lbl = cl->GetLabel(ilb); if (lbl<0 || lbl>=ntr) continue;
1692  int lblI = trIndArr[lbl];
1693  if (--lblI<0) continue; // not kept
1694  for (int icc=0;icc<kMaxCl;icc++) if (!clIndL[ilr][lblI+icc*ntrStore]) {clIndL[ilr][lblI+ntrStore*icc] = cl; break;} // first empty one
1695  }
1696  }
1697  }
1698  //
1699  Float_t clusterLay[2][AliITSMultReconstructor::kClNPar];
1700  double trComp[6][kMaxCl*kMaxCl];
1701  int indQual[kMaxCl*kMaxCl];
1702  //
1703  for (int itr=ntr;itr--;) {
1704  int lblI = trIndArr[itr];
1705  float weight = trWeightArr[itr];
1706  if (--lblI<0) continue; // discarded
1707  int ntrCand = 0; // number of tracklet candidates (including overlaps)
1708  for (int icl0=0;icl0<kMaxCl;icl0++) {
1709  AliITSRecPoint *cl0 = clIndL[0][lblI+icl0*ntrStore];
1710  if (!cl0 || !clIndL[1][lblI]) break;
1711  cl0->GetGlobalXYZ( clusterLay[0] );
1712  fMultReco->ClusterPos2Angles(clusterLay[0], fESDVtx);
1713  for (int icl1=0;icl1<kMaxCl;icl1++) {
1714  AliITSRecPoint *cl1 = clIndL[1][lblI+icl1*ntrStore];
1715  if (!cl1) break;
1716  cl1->GetGlobalXYZ( clusterLay[1] );
1717  fMultReco->ClusterPos2Angles(clusterLay[1], fESDVtx);
1718  trComp[AliITSMultReconstructor::kTrPhi][ntrCand] = clusterLay[0][AliITSMultReconstructor::kClPh];
1719  trComp[AliITSMultReconstructor::kTrTheta][ntrCand] = clusterLay[0][AliITSMultReconstructor::kClTh];
1720  trComp[AliITSMultReconstructor::kTrDTheta][ntrCand] = clusterLay[0][AliITSMultReconstructor::kClTh] - clusterLay[1][AliITSMultReconstructor::kClTh];
1721  trComp[AliITSMultReconstructor::kTrDPhi][ntrCand] = clusterLay[0][AliITSMultReconstructor::kClPh] - clusterLay[1][AliITSMultReconstructor::kClPh];
1722  trComp[AliITSMultReconstructor::kTrLab1][ntrCand] = icl1*10 + icl0;
1723  double &dphi = trComp[ntrCand][3];
1724  if (dphi>TMath::Pi()) dphi=2.*TMath::Pi()-dphi; // take into account boundary condition
1725  trComp[5][ntrCand] = fMultReco->CalcDist(trComp[AliITSMultReconstructor::kTrDPhi][ntrCand],
1726  trComp[AliITSMultReconstructor::kTrDTheta][ntrCand],
1727  trComp[AliITSMultReconstructor::kTrTheta][ntrCand]);
1728  ntrCand++;
1729  }
1730  }
1731  if (!ntrCand) continue; // no tracklets
1732  if (ntrCand>1) TMath::Sort(ntrCand,trComp[5],indQual,kFALSE); else indQual[0] = 0; // sort in weighted distance
1733  if (fRemoveOverlaps) ntrCand = 1; // select the best
1734  //
1735  // disable worst tracklet with shared cluster
1736  for (int itc=0;itc<ntrCand;itc++) {
1737  int ind = indQual[itc];
1738  if (trComp[AliITSMultReconstructor::kTrLab1][ind]<0) continue; // already disabled
1739  for (int jtc=itc+1;jtc<ntrCand;jtc++) {
1740  int jnd = indQual[jtc];
1741  if (trComp[AliITSMultReconstructor::kTrLab1][jnd]<0) continue; // already disabled
1742  if ( int(trComp[AliITSMultReconstructor::kTrLab1][ind])/10 == int(trComp[AliITSMultReconstructor::kTrLab1][jnd])/10 ||
1743  int(trComp[AliITSMultReconstructor::kTrLab1][ind])%10 == int(trComp[AliITSMultReconstructor::kTrLab1][jnd])%10) trComp[AliITSMultReconstructor::kTrLab1][jnd] = -1;
1744  }
1745  }
1746  //
1747  // store, but forbid cluster reusing
1748  TObjArray* histos = isPrimArr.TestBitNumber(itr) ? fHistosTrRcblPrim : fHistosTrRcblSec;
1749  for (int itc=0;itc<ntrCand;itc++) {
1750  int ind = indQual[itc];
1751  if (trComp[4][ind]<0) continue; // discarded
1752  double eta = -TMath::Log(TMath::Tan(trComp[AliITSMultReconstructor::kTrTheta][ind]/2));
1753  double phi = trComp[AliITSMultReconstructor::kTrPhi][ind]; // R+HACK
1754  if (eta<fEtaMin || eta>fEtaMax) continue;
1755  if (phi<fPhiMin || phi>fPhiMax) continue; // R+HACK
1756  double dThetaX = trComp[AliITSMultReconstructor::kTrTheta][ind];
1757  if (fScaleDTBySin2T) {
1758  double sint = TMath::Sin(trComp[AliITSMultReconstructor::kTrTheta][ind]);
1759  dThetaX /= (sint*sint);
1760  }
1761  FillHistosSet(histos,eta,
1762  //trComp[AliITSMultReconstructor::kTrPhi][ind],trComp[AliITSMultReconstructor::kTrTheta][ind],
1763  trComp[AliITSMultReconstructor::kTrDPhi][ind],trComp[AliITSMultReconstructor::kTrDTheta][ind],
1764  dThetaX,trComp[5][ind],weight);
1765  }
1766  }
1767  //
1768  delete[] clIndL[0];
1769  delete[] clIndL[1];
1770 }
1771 
1772 //_________________________________________________________________________
1774 {
1775  // fill info on clusters associated to good tracklets
1776  if (!fMultReco) return;
1777  int ntr = fMultReco->GetNTracklets();
1778  int clID[2];
1779  TH2F *hclU[2] = {(TH2F*)fHistosCustom->UncheckedAt(kHClUsedInfoL0),(TH2F*)fHistosCustom->UncheckedAt(kHClUsedInfoL1)};
1780  TH2F *hclA[2] = {(TH2F*)fHistosCustom->UncheckedAt(kHClAllInfoL0),(TH2F*)fHistosCustom->UncheckedAt(kHClAllInfoL1)};
1781  for (int itr=ntr;itr--;) {
1782  Float_t *trc = fMultReco->GetTracklet(itr);
1783  if (TMath::Abs(TMath::Abs(trc[AliITSMultReconstructor::kTrDPhi])-fDPhiShift)>fDPhiSCut) continue;
1784  if (fMultReco->CalcDist(trc[AliITSMultReconstructor::kTrDPhi],
1785  trc[AliITSMultReconstructor::kTrDTheta],
1786  trc[AliITSMultReconstructor::kTrTheta]) > fNStdCut) continue;
1787  clID[0] = (int)trc[AliITSMultReconstructor::kClID1];
1788  clID[1] = (int)trc[AliITSMultReconstructor::kClID2];
1789  for (int il=0;il<2;il++) {
1790  Float_t *clinf = fMultReco->GetClusterOfLayer(il,clID[il]);
1791  hclU[il]->Fill( clinf[AliITSMultReconstructor::kClZ], clinf[AliITSMultReconstructor::kClPh]);
1792  }
1793  }
1794  //
1795  for (int il=0;il<2;il++) for (int ic=fMultReco->GetNClustersLayer(il);ic--;) {
1796  Float_t *clinf = fMultReco->GetClusterOfLayer(il,ic);
1797  hclA[il]->Fill( clinf[AliITSMultReconstructor::kClZ], clinf[AliITSMultReconstructor::kClPh]);
1798  }
1799  //
1800 }
1801 
1802 //_________________________________________________________________________
1803 void AliTrackletTaskMulti::FillClusterInfoFromMult(const AliMultiplicity* mlt, double zVertex)
1804 {
1805  // fill info on clusters taking them from Multiplicity object
1806  const double kRSPD2 = 3.9;
1807  TH2F *hclU = (TH2F*)fHistosCustom->UncheckedAt(kHClUsedInfoL0);
1808  TH2F *hclA = (TH2F*)fHistosCustom->UncheckedAt(kHClAllInfoL0);
1809  int ntr = mlt->GetNumberOfTracklets();
1810  for (int itr=ntr;itr--;) {
1811  Bool_t goodTracklet = kTRUE;
1812  if (TMath::Abs( mlt->GetDeltaPhi(itr)-fDPhiShift)>fDPhiSCut) goodTracklet = kFALSE;
1813  if (mlt->CalcDist(itr) > fNStdCut) goodTracklet = kFALSE;
1814  double phi = mlt->GetPhi(itr);
1815  double z = kRSPD2/TMath::Tan(mlt->GetTheta(itr)) + zVertex;
1816  if (goodTracklet) hclU->Fill(z,phi);
1817  hclA->Fill(z,phi);
1818  }
1819  //
1820  int ncl = mlt->GetNumberOfSingleClusters();
1821  for (int icl=ncl;icl--;) {
1822  double phi = mlt->GetPhiSingle(icl);
1823  double z = kRSPD2/TMath::Tan(mlt->GetThetaSingle(icl)) + zVertex;
1824  hclA->Fill(z,phi);
1825  }
1826  //
1827 }
1828 
1829 //______________________________________________
1831 {
1832  int nv = sizeof(fgCentSelName)/sizeof(char*);
1833  for (int i=0;i<nv;i++) {
1834  if (!strcmp(var, fgCentSelName[i])) {
1835  return;
1836  }
1837  }
1838  AliFatalF("Unknown centrality var: %s",var);
1839 }
1840 
1841 //______________________________________________
1843 {
1844  // set user defined percentiles
1845  fCentPerc.Set(nbins+1);
1846  fNCentBins = nbins;
1847  AliInfoF("Defining %d centrality bins",fNCentBins);
1848  for (int i=0;i<=nbins;i++) {
1849  fCentPerc[i] = arr[i];
1850  if (i<nbins) AliInfoF("CentBin# %.1f-%.1f",arr[i],arr[i+1]);
1851  }
1852  //
1853 }
1854 
1855 //______________________________________________
1857 {
1858  // set user defined percentiles
1859  fCentPerc.Set(nbins+1);
1860  fNCentBins = nbins;
1861  AliInfoF("Defining %d centrality bins",fNCentBins);
1862  for (int i=0;i<=nbins;i++) {
1863  fCentPerc[i] = arr[i];
1864  if (i<nbins) AliInfoF("CentBin# %.1f-%.1f",arr[i],arr[i+1]);
1865  }
1866  //
1867 }
1868 
1869 //______________________________________________
1870 // R+HACK
1871 //
1872 TFile *file_REWEIGHTpt = NULL;
1874 const Char_t *h_REWEIGHTpt_name[11] = {"c0_5", "c0_5", "c5_10", "c5_10", "c10_20", "c20_30", "c30_40", "c40_50", "c50_60", "c60_70", "c70_80"};
1875 //
1876 TFile *file_REWEIGHTpid = NULL;
1880 //
1881 TFile *file_REWEIGHTstr = NULL;
1889 //
1891 {
1892  // reweight particle weights in MC stack
1893  if (!fStack)
1894  return;
1895 
1896  if (!file_REWEIGHTpt)
1897  file_REWEIGHTpt = TFile::Open("REWEIGHTpt.root");
1898  h_REWEIGHTpt = 0;
1899  Double_t cc = ((TH1*)fHistosCustom->UncheckedAt(kHStatCent))
1900  ->GetXaxis()->GetBinCenter(fCurrCentBin+1);
1901  Int_t c1 = 0, c2 = 5;
1902  if (cc <= 5) { c1 = 0; c2 = 5; }
1903  else if (cc <= 10) { c1 = 5; c2 = 10; }
1904  else if (cc <= 20) { c1 = 10; c2 = 20; }
1905  else if (cc <= 30) { c1 = 20; c2 = 30; }
1906  else if (cc <= 40) { c1 = 30; c2 = 40; }
1907  else if (cc <= 50) { c1 = 40; c2 = 50; }
1908  else if (cc <= 60) { c1 = 50; c2 = 60; }
1909  else if (cc <= 70) { c1 = 60; c2 = 70; }
1910  else { c1 = 70; c2 = 80; }
1911  TString ptName; ptName.Form("ptWeight_c%d_%d", c1, c2);
1912  h_REWEIGHTpt = (TH1 *)file_REWEIGHTpt->Get(ptName);
1913  if (!h_REWEIGHTpt) {
1914  Warning("ReweightStack", "pT weight histogram %s (%d) not found in %s",
1915  ptName.Data(), fCurrCentBin, file_REWEIGHTpt->GetName());
1916  }
1917 
1918  if (!file_REWEIGHTpid) {
1919  switch (fReweightFlag) {
1920  case 0:
1921  file_REWEIGHTpid = TFile::Open("REWEIGHTpid.root");
1922  break;
1923  case 1:
1924  file_REWEIGHTpid = TFile::Open("REWEIGHTpid_pi+.root");
1925  break;
1926  case -1:
1927  file_REWEIGHTpid = TFile::Open("REWEIGHTpid_pi-.root");
1928  break;
1929  case 2:
1930  file_REWEIGHTpid = TFile::Open("REWEIGHTpid_ka+.root");
1931  break;
1932  case -2:
1933  file_REWEIGHTpid = TFile::Open("REWEIGHTpid_ka-.root");
1934  break;
1935  case 3:
1936  file_REWEIGHTpid = TFile::Open("REWEIGHTpid_pr+.root");
1937  break;
1938  case -3:
1939  file_REWEIGHTpid = TFile::Open("REWEIGHTpid_pr-.root");
1940  break;
1941  }
1942  h_REWEIGHTpid_pi = (TH1 *)file_REWEIGHTpid->Get("pidWeight_pi");
1943  h_REWEIGHTpid_ka = (TH1 *)file_REWEIGHTpid->Get("pidWeight_ka");
1944  h_REWEIGHTpid_pr = (TH1 *)file_REWEIGHTpid->Get("pidWeight_pr");
1945  }
1946 
1947  if (!file_REWEIGHTstr) {
1948  switch (fReweightFlag) {
1949  case 0:
1950  file_REWEIGHTstr = TFile::Open("REWEIGHTstr.root");
1951  break;
1952  case 1:
1953  file_REWEIGHTstr = TFile::Open("REWEIGHTstr+.root");
1954  break;
1955  case -1:
1956  file_REWEIGHTstr = TFile::Open("REWEIGHTstr-.root");
1957  break;
1958  }
1959  h_REWEIGHTstr_pi = (TH1 *)file_REWEIGHTstr->Get("strWeight_pi");
1960  h_REWEIGHTstr_ka = (TH1 *)file_REWEIGHTstr->Get("strWeight_ka");
1961  h_REWEIGHTstr_pr = (TH1 *)file_REWEIGHTstr->Get("strWeight_pr");
1962  h_REWEIGHTstr_k0 = (TH1 *)file_REWEIGHTstr->Get("strWeight_k0");
1963  h_REWEIGHTstr_la = (TH1 *)file_REWEIGHTstr->Get("strWeight_la");
1964  h_REWEIGHTstr_si = (TH1 *)file_REWEIGHTstr->Get("strWeight_si");
1965  h_REWEIGHTstr_xi = (TH1 *)file_REWEIGHTstr->Get("strWeight_xi");
1966  }
1967 
1968  int ntr = fStack->GetNtrack();
1969  for (int itr=ntr;itr--;) {
1970  TParticle* part = fStack->Particle(itr);
1971  if (!part) continue;
1972  part->SetWeight(GetPrimaryWeight(itr));
1973  }
1974 }
1975 
1976 //______________________________________________
1977 // R+HACK
1979 {
1980  // get weight for a primary particle
1981  // if secondary, go back till primary mother
1982 
1983  if (!fStack || lab < 0)
1984  return 1.;
1985 
1986  TParticle *part = fStack->Particle(lab);
1987  // if not primary, iterate on first mother
1988  if (!fStack->IsPhysicalPrimary(lab)) {
1989  int lmom = part->GetFirstMother();
1990  return GetPrimaryWeight(lmom);
1991  }
1992 
1993  // this is a primary particle
1994  float weight = part->GetWeight();
1995  float pt = part->Pt();
1996  int pdg = TMath::Abs(part->GetPdgCode());
1997  if ((fReweightStack & kReweightPt) && pt < 5.) {
1998  if (h_REWEIGHTpt) {
1999  weight *= h_REWEIGHTpt->GetBinContent(h_REWEIGHTpt->FindBin(pt));
2000  if (pt < 0.05) {
2001  switch (fReweightFlag) {
2002  case 1:
2003  weight *= 1.3;
2004  break;
2005  case -1:
2006  weight *= 0.7;
2007  break;
2008  }
2009  }
2010  // printf(">>> REWEIGHTpt: pt=%f, w=%f\n", pt, weight);
2011  }
2012  }
2013  if ((fReweightStack & kReweightPid)) {
2014  switch (pdg) {
2015  case 211:
2016  weight *= h_REWEIGHTpid_pi->GetBinContent(1);
2017  break;
2018  case 2212:
2019  weight *= h_REWEIGHTpid_pr->GetBinContent(1);
2020  break;
2021  case 321:
2022  weight *= h_REWEIGHTpid_ka->GetBinContent(1);
2023  break;
2024  }
2025  // printf(">>> REWEIGHTpid: pdg=%d, w=%f\n", pdg, weight);
2026  }
2027  if ((fReweightStack & kReweightStr)) {
2028  switch (pdg) {
2029  case 211:
2030  weight *= h_REWEIGHTstr_pi->GetBinContent(1);
2031  break;
2032  case 2212:
2033  weight *= h_REWEIGHTstr_pr->GetBinContent(1);
2034  break;
2035  case 321:
2036  weight *= h_REWEIGHTstr_ka->GetBinContent(1);
2037  break;
2038  case 310:
2039  weight *= h_REWEIGHTstr_k0->GetBinContent(1);
2040  break;
2041  case 3122:
2042  weight *= h_REWEIGHTstr_la->GetBinContent(1);
2043  break;
2044  case 3112: case 3222:
2045  weight *= h_REWEIGHTstr_si->GetBinContent(1);
2046  break;
2047  case 3312:
2048  weight *= h_REWEIGHTstr_xi->GetBinContent(1);
2049  break;
2050  }
2051  // printf(">>> REWEIGHTpid: pdg=%d, w=%f\n", pdg, weight);
2052  }
2053 
2054  return weight;
2055 }
TObjArray * fHistosTrInj
all tracklets in data
Int_t pdg
const Color_t cc[]
Definition: DrawKs.C:1
TH1 * h_REWEIGHTpid_pi
Bool_t HaveCommonParent(const float *clLabs0, const float *clLabs1)
double Double_t
Definition: External.C:58
Definition: External.C:260
TH1 * h_REWEIGHTstr_la
TFile * file_REWEIGHTpt
Definition: External.C:236
Float_t fNPart
N of signal tracklets |eta|<2.
Int_t GetPdgBin(Int_t pdgCode)
void SetRemoveOverlaps(Bool_t w=kFALSE)
Bool_t FillHistosSet(TObjArray *histos, double eta, double dphi, double dtheta, double dthetaX, double dist, double weight)
static const char * fgCentSelName[]
static const char * fgkPDGNames[]
centrality types
void Print(Option_t *option) const
Definition: External.C:244
void AddHisto(TObjArray *histos, TObject *h, Int_t at=-1)
Int_t GetCentralityBin(Float_t percentile) const
TH1 * h_REWEIGHTstr_k0
void SetThetaWindow(float w=0.025)
TH1 * h_REWEIGHTpt
TObjArray * fHistosTrSec
primary
char Char_t
Definition: External.C:18
Float_t fVtxMC[3]
ESD vertex.
TObjArray * fHistosCustom
Secondary Reconstructable.
void SetPhiWindow(float w=0.08)
Bool_t fUseSpecialOutput
rec.vertex is good
Bool_t GetDoInjection() const
TH1 * h_REWEIGHTstr_si
Int_t fNTreta2
N of primaries |eta|<2.
TString kData
Declare data MC or deltaAOD.
Bool_t GetDoNormalReco() const
TH1 * h_REWEIGHTstr_ka
AliTrackletTaskMulti(const char *name="AliTrackletTaskMulti")
TH1 * h_REWEIGHTpid_pr
AliITSMultRecBg * fMultReco
TObjArray * fHistosTrPrim
rotated
int Int_t
Definition: External.C:63
void CheckCentralityVar(const char *var)
unsigned int UInt_t
Definition: External.C:33
float Float_t
Definition: External.C:68
const double kEtaBinWidth
void SetCentPercentiles(Float_t *arr, Int_t nbins)
TH1 * h_REWEIGHTstr_pr
TFile * file_REWEIGHTstr
TObjArray * fHistosTrRot
injected
virtual void UserExec(Option_t *option)
AliMCEvent * fMCEvent
MC stack.
TObjArray * fHistosTrCombU
combinatorials
static const Int_t fgkPDGCodes[]
pdg names
AliStack * fStack
tree of recpoints
void FillHistos(Int_t type, const AliMultiplicity *mlt)
void FillSpecies(Int_t primsec, Int_t id)
TH1 * h_REWEIGHTstr_pi
Bool_t fVtxOK
did current event pass phys.sel.?
Float_t fEtaMin
custom histos
Int_t fNPrimMCeta2
MC gen vertex.
void ReweightStack(double weight)
TTree * fRPTree
mult.reco object
TObjArray * BookHistosSet(const char *pref, UInt_t selHistos=0xffffffff)
Definition: External.C:220
const char Option_t
Definition: External.C:48
void FillClusterInfoFromMult(const AliMultiplicity *mlt, double zVertex)
void SetNStdDev(Float_t f=1.)
void SetPhiOverlapCut(float w=0.005)
Float_t fESDVtx[3]
MC Event.
TObjArray * fHistosTrRcblSec
Primary Reconstructable.
void SetPhiShift(float w=0.0045)
const Int_t nbins
TObjArray * fHistosTrRcblPrim
combinatorials uncorrelated
bool Bool_t
Definition: External.C:53
TFile * file_REWEIGHTpid
virtual void Terminate(Option_t *)
TH1 * h_REWEIGHTstr_xi
TObjArray * fHistosTrComb
secondary
virtual void UserCreateOutputObjects()
Bool_t GetDoRotation() const
void SetZetaOverlapCut(float w=0.05)
TH1 * h_REWEIGHTpid_ka
void SetScaleDThetaBySin2T(Bool_t v=kFALSE)
const Char_t * h_REWEIGHTpt_name[11]
Definition: External.C:196
TList * OpenFile(const char *fname)
Definition: DrawAnaELoss.C:65