]> git.uio.no Git - u/mrichter/AliRoot.git/blame - PWGUD/multVScentPbPb/AliTrackletTaskUni.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Fix typo
[u/mrichter/AliRoot.git] / PWGUD / multVScentPbPb / AliTrackletTaskUni.cxx
CommitLineData
a9a39f46 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"
3c78f321 42//#include "../ANALYSIS/EventMixing/AliMixEventInputHandler.h"
a9a39f46 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"
c79a7c1c 60#include "../ITS/AliITSsegmentationSPD.h"
a9a39f46 61#include "AliLog.h"
62
63#include "AliPhysicsSelection.h"
a9a39f46 64#include "AliTrackletTaskUni.h"
65#include "AliITSMultRecBg.h"
66#include "AliGenEventHeader.h"
67#include "AliGenHijingEventHeader.h"
68#include "AliGenDPMjetEventHeader.h"
69
70ClassImp(AliTrackletTaskUni)
71
72const 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
123const 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
176AliTrackletTaskUni::AliTrackletTaskUni(const char *name)
177 : AliAnalysisTaskSE(name),
178*/
179//________________________________________________________________________
180AliTrackletTaskUni::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//
3c78f321 207 fEtaMin(-3.0),
208 fEtaMax(3.0),
a9a39f46 209 fZVertexMin(-20),
210 fZVertexMax( 20),
211 fMultCutMin(0),
212 fMultCutMax(99999),
3c78f321 213 fMCV0Scale(0.7520),
a9a39f46 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),
c79a7c1c 231 fMCevent(0),
3c78f321 232 fDontMerge(kFALSE)
a9a39f46 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//________________________________________________________________________
258AliTrackletTaskUni::~AliTrackletTaskUni()
259{
260 // Destructor
261 // histograms are in the output list and deleted when the output
262 // list is deleted by the TSelector dtor
c79a7c1c 263 if (!AliAnalysisManager::GetAnalysisManager()->IsProofMode()) { //RRR
a9a39f46 264 printf("Deleteing output\n");
265 delete fOutput;
a9a39f46 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//________________________________________________________________________
285void AliTrackletTaskUni::UserCreateOutputObjects()
286{
c79a7c1c 287 // create output
3c78f321 288
289 if (fDontMerge) {
290 OpenFile(1);
291 printf("No merging will be done\n");
292 }
293
a9a39f46 294 fOutput = new TList();
295 fOutput->SetOwner();
296 //
297 AliCDBManager *man = AliCDBManager::Instance();
298 if (fUseMC) {
3c78f321 299 printf("Loading MC geometry\n");
a9a39f46 300 Bool_t newGeom = kTRUE;
301 man->SetDefaultStorage("alien://Folder=/alice/simulation/2008/v4-15-Release/Residual");
302 if (newGeom) {
303 // new geom
3c78f321 304 AliCDBEntry* obj = man->Get("GRP/Geometry/Data",130844,-1,-1);
a9a39f46 305 AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
3c78f321 306 if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",130844,-1,-1)) AliFatal("Failed to misalign geometry");
a9a39f46 307 }
308 else {
309 // old geom
3c78f321 310 AliCDBEntry* obj = man->Get("GRP/Geometry/Data",130844,7,-1);
a9a39f46 311 AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
312 if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",130845,5,-1)) AliFatal("Failed to misalign geometry");
313 }
314 }
315 else {
3c78f321 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);
a9a39f46 319 AliGeomManager::SetGeometry((TGeoManager*) obj->GetObject());
3c78f321 320 if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",137366,-1,-1)) AliFatal("Failed to misalign geometry");
a9a39f46 321 }
322 //
3c78f321 323 printf("Geometry loaded\n");
a9a39f46 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
3c78f321 371//________________________________________________________________________
372void AliTrackletTaskUni::RegisterStat()
373{
c79a7c1c 374 // account conditions
3c78f321 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
a9a39f46 410//________________________________________________________________________
411void AliTrackletTaskUni::UserExec(Option_t *)
412{
413 // Main loop
414 //
3c78f321 415 if (fDontMerge) RegisterStat();
416 Bool_t needRecPoints = GetDoNormalReco() || GetDoInjection() || GetDoRotation() || GetDoMixing();
417
a9a39f46 418 AliAnalysisManager* anMan = AliAnalysisManager::GetAnalysisManager();
419 fRPTree = fRPTreeMix = 0;
3c78f321 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 }
a9a39f46 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 */
3c78f321 452 TH1F* hstat = (TH1F*)fHistosCustom->UncheckedAt(kHStat);
453 //
454 hstat->Fill(kEvTot0); // RS
a9a39f46 455 const AliESDVertex* vtxESD = esd->GetPrimaryVertexSPD();
456 if (vtxESD->GetNContributors()<1) return;
3c78f321 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 //
a9a39f46 468 const AliMultiplicity* multESD = esd->GetMultiplicity();
3c78f321 469 /*
a9a39f46 470 const AliESDVertex* vtxESDTPC = esd->GetPrimaryVertexTPC();
471 if (vtxESDTPC->GetNContributors()<1 ||
472 vtxESDTPC->GetNContributors()<(-10.+0.25*multESD->GetNumberOfITSClusters(0))) return;
3c78f321 473 */
a9a39f46 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
3c78f321 485 //------------------------------------------------------
a9a39f46 486 AliESDZDC *esdZDC = esd->GetESDZDC();
3c78f321 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.;
a9a39f46 491 }
a9a39f46 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) {
3c78f321 500 const double v0Scale = fMCV0Scale;
a9a39f46 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 //
3c78f321 525 // printf("Multiplicity from ESD:\n");
a9a39f46 526 // multESD->Print();
527 //
528 AliMCEventHandler* eventHandler = 0;
c79a7c1c 529 fMCevent = 0;
a9a39f46 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; }
c79a7c1c 535 fMCevent = eventHandler->MCEvent();
536 if (!fMCevent) { printf("ERROR: Could not retrieve MC event\n"); return; }
537 fStack = fMCevent->Stack();
a9a39f46 538 if (!fStack) { printf("Stack not available\n"); return; }
539 }
540 //
3c78f321 541 if (needRecPoints) {
542 fRPTree = handRP->GetTreeR("ITS");
543 if (!fRPTree) { AliError(" Invalid ITS cluster tree !\n"); return; }
544 }
a9a39f46 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
3c78f321 558 static int prnRec = 10;
559 static int prnInj = 10;
560 //
a9a39f46 561 hstat->Fill(kEvProcData);
562 if (GetDoNormalReco() || GetDoInjection()) { // for the injection the normal reco should be done
563 InitMultReco();
564 fMultReco->Run(fRPTree, vtxf);
a9a39f46 565 AliMultiplicity* mlt = fMultReco->GetMultiplicity();
3c78f321 566 if (mlt && (--prnRec)>0) {
567 printf("Multiplicity Reconstructed: %d\n",prnRec);
568 mlt->Print();
569 }
a9a39f46 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);
3c78f321 579 fMultReco->Run(fRPTree, vtxf);
a9a39f46 580 AliMultiplicity* mlt = fMultReco->GetMultiplicity();
3c78f321 581 if (mlt && (--prnInj)>0) {
582 printf("Multiplicity from Injection: %d\n",prnInj);
583 mlt->Print();
584 }
585 // if (mlt) mlt->Print();
a9a39f46 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);
a9a39f46 596 AliMultiplicity* mlt = fMultReco->GetMultiplicity();
a9a39f46 597 hstat->Fill(kEvProcRot);
598 FillHistos(kBgRot,mlt);
599 }
600 //
601 if (GetDoMixing()) {
3c78f321 602 /*
a9a39f46 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);
3c78f321 627 */
628 AliFatal("Mixing is outphased");
a9a39f46 629 //
630 }
631 // =============================================================================<<<
632 //
633 delete fMultReco;
634 fMultReco = 0;
635 //
636}
637//________________________________________________________________________
638void AliTrackletTaskUni::Terminate(Option_t *)
639{
c79a7c1c 640 // print itself
a9a39f46 641 Printf("Terminating...");
3c78f321 642 RegisterStat();
a9a39f46 643 // AliAnalysisTaskSE::Terminate();
644}
645
646
647//_________________________________________________________________________
648void 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//_________________________________________________________________________
667TObjArray* 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 ------------------------------>>>
3c78f321 674 hstat = new TH1F("hStat","Run statistics",kNStatBins,0.5,kNStatBins+0.5);
a9a39f46 675 //
3c78f321 676 for (int ib=1;ib<=kNStatBins;ib++) hstat->GetXaxis()->SetBinLabel(ib,"-"); // dummy label
677 hstat->GetXaxis()->SetBinLabel(kEvTot0, "Ev.Tot0");
a9a39f46 678 hstat->GetXaxis()->SetBinLabel(kEvTot, "Ev.Tot");
3c78f321 679 hstat->GetXaxis()->SetBinLabel(kEvTotPlp, "Ev.Tot Plp");
680 //
a9a39f46 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");
3c78f321 697 hstat->GetXaxis()->SetBinLabel(kEtaMin,"#eta_{min}");
698 hstat->GetXaxis()->SetBinLabel(kEtaMax,"#eta_{max}");
a9a39f46 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 //----------------------------------------------------------------------
3c78f321 782 int nEtaBinsS = int((fEtaMax-fEtaMin)/0.1);
a9a39f46 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
3c78f321 790 TH2F* hzvetap = new TH2F("zvEtaPrimMC","Z vertex vs eta PrimMC",nEtaBinsS,fEtaMin,fEtaMax,nZVBinsS,fZVertexMin,fZVertexMax);
a9a39f46 791 hzvetap->GetXaxis()->SetTitle("#eta");
792 hzvetap->GetYaxis()->SetTitle("Zvertex");
793 AddHisto(histos,hzvetap,kHZVEtaPrimMC);
794 }
795 //
796 if (GetDoMixing()) {
797 //
3c78f321 798 /*
a9a39f46 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);
3c78f321 809 */
a9a39f46 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//_________________________________________________________________________
837TObjArray* 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;
3c78f321 843 int nDistBins = int(fNStdDev)*4;
a9a39f46 844
3c78f321 845 int nEtaBins = int((fEtaMax-fEtaMin)/0.1);
a9a39f46 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};
3c78f321 863 double xmnEZD[3] = { fEtaMin,fZVertexMin,0};
864 double xmxEZD[3] = { fEtaMax,fZVertexMax,fNStdDev};
a9a39f46 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)};
3c78f321 879 double xmnEZP[3] = { fEtaMin,fZVertexMin,-dphir};
880 double xmxEZP[3] = { fEtaMax,fZVertexMax, dphir};
a9a39f46 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);
3c78f321 891 h2 = new TH2F(buffn,bufft,nEtaBins,fEtaMin,fEtaMax, nZVBins, fZVertexMin,fZVertexMax);
a9a39f46 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);
3c78f321 919 h2 = new TH2F(buffn,bufft,nEtaBins, fEtaMin,fEtaMax,kNDPhiBins,-dphir,dphir);
a9a39f46 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)":"");
3c78f321 928 h2 = new TH2F(buffn,bufft,nEtaBins, fEtaMin,fEtaMax,kNDThtBins,-dthtr,dthtr);
a9a39f46 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);
3c78f321 938 h2 = new TH2F(buffn,bufft,nEtaBins, fEtaMin,fEtaMax,nDistBins,0,fNStdDev);
a9a39f46 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//_________________________________________________________________________
980void 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//_________________________________________________________________________
989void 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) {
c79a7c1c 1007 pyHeader = (AliGenHijingEventHeader*) fMCevent->GenEventHeader();//header->GenEventHeader();
a9a39f46 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//_________________________________________________________________________
1072void AliTrackletTaskUni::FillMCPrimaries(TH2F* hetaz)
1073{
1074 // fill all MC primaries Zv vs Eta
c79a7c1c 1075 if (!fStack || !fMCevent) return;
a9a39f46 1076
1077 //---------------------------------------- CHECK ------------------------------>>>
1078 TArrayF vtxMC;
1079 AliGenHijingEventHeader* pyHeader = 0;
1080 //
1081 if (fUseMC) {
c79a7c1c 1082 pyHeader = (AliGenHijingEventHeader*) fMCevent->GenEventHeader();//header->GenEventHeader();
a9a39f46 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;
c79a7c1c 1092 AliMCParticle *part = (AliMCParticle*)fMCevent->GetTrack(itr);
a9a39f46 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//_________________________________________________________________________
1110void 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));
3c78f321 1119 if (eta<fEtaMin || eta>fEtaMax) return;
a9a39f46 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 //
3c78f321 1139 if (histos->UncheckedAt(kHDPhiDTheta))
1140 ((TH2F*)histos->UncheckedAt(kHDPhiDTheta))->Fill(dphi,dtheta);
1141 //
a9a39f46 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//__________________________________________________________________
1169void 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//_________________________________________________________________________
1202Int_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//_________________________________________________________________________
1216Bool_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();
3c78f321 1224 // printf("\nHaveCommonParent \n");
a9a39f46 1225 for (int il=0;il<2;il++) {
3c78f321 1226
a9a39f46 1227 for (int ilb=0;ilb<3;ilb++) {
1228 int lbl = (int)labs[il][ilb];
3c78f321 1229 if (lbl<2 || lbl>=ntr) continue;
a9a39f46 1230 //
1231 while (npars[il]<kMaxPar-1) {
a9a39f46 1232 TParticle* part = fStack->Particle(lbl);
1233 if (!part) break;
3c78f321 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;
a9a39f46 1239 lbl = part->GetFirstMother();
1240 if (lbl<1 || lbl>=ntr) break;
1241 }
3c78f321 1242 // printf("\n");
a9a39f46 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//_________________________________________________________________________
1252void 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