AliPhysics  db95e02 (db95e02)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliTrackletAODTask.C
Go to the documentation of this file.
1 
11 #include <AliAnalysisTaskSE.h>
12 #include <TParameter.h>
13 #ifndef __CINT__
14 #include "AliAODTracklet.C"
15 #include "AliVVertex.h"
16 #include <AliVEvent.h>
17 #include <AliAODEvent.h>
18 #include <AliVVertex.h>
19 #include <AliVertex.h>
20 #include <AliVMultiplicity.h>
21 #include <AliMultiplicity.h>
22 #include <AliAnalysisManager.h>
23 #include <AliVEventHandler.h>
24 #include <AliESDInputHandlerRP.h>
25 #include <AliITSMultRecBg.h>
26 #include <AliCDBManager.h>
27 #include <AliCDBEntry.h>
28 #include <AliCDBPath.h>
29 #include <AliCDBId.h>
30 #include <AliGeomManager.h>
31 #include <AliAODHandler.h>
32 #include <TGeoManager.h>
33 #include <TGeoGlobalMagField.h>
34 #include <climits>
35 #else
36 class AliAODTracklet;
37 class AliVEvent;
38 class AliVVertex;
39 class AliMultiplicity;
40 class AliVMultiplicity;
41 class AliAnalysisDataContainer;
42 class AliITSMultRecBg;
43 class AliGeomManager; // Auto-load
44 class AliCDBManager; // Auto-load
45 class AliMultSelection; // Auto-load
46 class TClonesArray;
47 class TTree;
48 class TGeoGlobalMagField; // Auto-load
49 class TGeoManager; // Auto-load
50 class TParticle;
51 #endif
52 
53 
54 //====================================================================
61 {
62 public:
72  AliTrackletAODTask(const char* name);
90  void Print(Option_t*) const;
96  static AliTrackletAODTask* Create();
109  void UserExec(Option_t*);
114  void FinishTaskOutput() { /*WorkerFinalize();*/ }
119  void Terminate(Option_t*) {}
125  Bool_t Connect();
126  /* @} */
136  void SetScaleDTheta(Bool_t x=false) { fScaleDTheta = x; }
142  void SetDPhiShift(Double_t x=0.0045) { fDPhiShift = x; }
148  void SetMaxDelta(Double_t x=25) { fMaxDelta = x; }
154  void SetDThetaWindow(Double_t x=0.025) { fDThetaWindow = x; }
160  void SetDPhiWindow(Double_t x=0.06) { fDPhiWindow = x; }
166  void SetPhiOverlapCut(Double_t x=0.005) { fPhiOverlapCut = x; }
183  /* @} */
184 protected:
192  virtual Bool_t WorkerInit();
198  Bool_t InitCDB();
204  virtual Int_t GetCDBReferenceRun() const { return 137161; }
210  virtual const char* GetCDBReferenceURL() const
211  {
212  return "alien://Folder=/alice/data/2010/OCDB";
213  }
219  Bool_t InitBranch();
230  virtual Bool_t ProcessEvent();
242  Bool_t HasField(AliVEvent* event);
248  TTree* FindClusters();
254  virtual TTree* FilterClusters(TTree* t) { return t; }
260  virtual void CleanClusters(TTree*& t) {};
266  AliVEvent* FindEvent();
276  const AliVVertex* FindIP(AliVEvent* event,
277  Double_t maxDispersion=0.04,
278  Double_t maxZError=0.25);
285  Bool_t Reconstruct(TTree* clusters, const AliVVertex* ip);
292  /* @} */
302  virtual const char* TrackletClassName() const { return "AliAODTracklet"; }
310  virtual AliAODTracklet* MakeTracklet(Bool_t normal);
319  Bool_t ProcessTracklets(Bool_t normal, AliVMultiplicity* mult);
329  virtual AliAODTracklet* ProcessTracklet(Bool_t normal,
330  AliMultiplicity* mult,
331  Int_t no);
332  /* @} */
333 
335  TClonesArray* fTracklets;
336  /* Output container */
337  // Container* fContainer;
353  AliITSMultRecBg* fReco;
354 
358 
360 };
361 //____________________________________________________________________
363  : AliAnalysisTaskSE(),
364  fTracklets(0),
365  // fContainer(0),
366  fScaleDTheta(true),
367  fMaxDelta(25),
368  fDPhiShift(0.0045),
369  fDThetaWindow(0.025),
370  fDPhiWindow(0.06),
371  fPhiOverlapCut(0.005),
372  fZEtaOverlapCut(0.05),
373  fReco(0),
374  fFilterStrange(0),
375  fStrangeLoss(0)
376 {}
377 //____________________________________________________________________
379  : AliAnalysisTaskSE("tracklets"),
380  fTracklets(0),
381  // fContainer(0),
382  fScaleDTheta(true),
383  fMaxDelta(25),
384  fDPhiShift(0.0045),
385  fDThetaWindow(0.025),
386  fDPhiWindow(0.06),
387  fPhiOverlapCut(0.005),
388  fZEtaOverlapCut(0.05),
389  fReco(0),
390  fFilterStrange(0),
391  fStrangeLoss(0)
392 {
393  // DefineOutput(1,TList::Class());
394 }
395 //____________________________________________________________________
397  : AliAnalysisTaskSE(other),
398  fTracklets(0),
399  // fContainer(0),
400  fScaleDTheta(other.fScaleDTheta),
401  fMaxDelta(other.fMaxDelta),
402  fDPhiShift(other.fDPhiShift),
403  fDThetaWindow(other.fDThetaWindow),
404  fDPhiWindow(other.fDPhiWindow),
405  fPhiOverlapCut(other.fPhiOverlapCut),
406  fZEtaOverlapCut(other.fZEtaOverlapCut),
407  fReco(0),
408  fFilterStrange(other.fFilterStrange),
409  fStrangeLoss(0)
410 {}
411 //____________________________________________________________________
413 {
414  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
415  if (!mgr) {
416  AliError("No analysis manager to connect to.");
417  return false;
418  }
419  AliAODHandler* ah =
420  dynamic_cast<AliAODHandler*>(mgr->GetOutputEventHandler());
421  if (!ah) {
422  AliError("No AOD output handler!");
423  return false;
424  }
425 
426  // Add to the manager
427  mgr->AddTask(this);
428 
429  // Always connect input
430  mgr->ConnectInput(this, 0, mgr->GetCommonInputContainer());
431 
432  return true;
433 }
434 //____________________________________________________________________
436 {
437  Printf("%s: %s", ClassName(), GetName());
438  Printf(" %22s: %d", "Scale by sin^2(theta)", fScaleDTheta);
439  Printf(" %22s: %f", "Delta phi shift", fDPhiShift);
440  Printf(" %22s: %f", "max Delta", fMaxDelta);
441  Printf(" %22s: %f", "Delta theta window", fDThetaWindow);
442  Printf(" %22s: %f", "Delta phi window", fDPhiWindow);
443  Printf(" %22s: %f", "phi overlap cut", fPhiOverlapCut);
444  Printf(" %22s: %f", "z-eta overlap cut", fZEtaOverlapCut);
445  Printf(" %22s: %s", "Filter strange",
446  fFilterStrange==0 ? "no" :
447  fFilterStrange==1 ? "random" : "track");
448 }
449 
450 //____________________________________________________________________
452 {
453  return *this;
454 }
455 //____________________________________________________________________
457 {
458  // fContainer = new Container;
459  // fContainer->SetOwner();
460  // fContainer->SetName(GetName());
461  // PostData(1, fContainer);
462 
463  // Create container of tracklets
464  if (WorkerInit()) return;
465 
466  AliWarning("Failed to initialize task on worker, making zombie");
467  SetZombie();
468 
469 }
470 
471 //____________________________________________________________________
473 {
474  AliInfo("Now initialising CDB");
475  AliAnalysisManager* anaMgr = AliAnalysisManager::GetAnalysisManager();
476  if (!anaMgr) {
477  AliError("No manager defined!");
478  return false;
479  }
480  // Check if we have the CDB connect task, and if so, do nothing as
481  // we rely on that task set up things properly.
482  const char* cdbNames[] = {"CDBconnect", "cdb", 0 };
483  const char** ptr = cdbNames;
484  while (*ptr) {
485  AliAnalysisTask* cdbConnect = anaMgr->GetTask(*ptr);
486  if (cdbConnect && cdbConnect->IsA()->InheritsFrom("AliTaskCDBconnect")) {
487  AliInfoF("CDB-connect task (%s: %s) present, do nothing",
488  cdbConnect->ClassName(), *ptr);
489  return true;
490  }
491  ptr++;
492  }
493  // Otherwise, we need to do stuff ourselves
494  AliInfo("Get the CDB manager");
495  AliCDBManager* cdbMgr = AliCDBManager::Instance();
496  if (!cdbMgr) {
497  AliError("Failed to get instance of CDB manager");
498  return false;
499  }
500  Int_t refRun = GetCDBReferenceRun();
501  TString refUrl = GetCDBReferenceURL();
502  AliWarningF("Using reference CDB storage \"%s\" and run \"%d\"",
503  refUrl.Data(), refRun);
504  // Set a very particular default storage. Perhaps we can do this
505  // with specific storages instead! Depends on whether the
506  // reconstruction also uses CDB - probably does - in which case we
507  // need specific storages for that too.
508  cdbMgr->SetDefaultStorage(refUrl);
509  // Now load our geometry - from a LHC10h run
510  AliInfo("Get Geometry entry");
511  AliCDBEntry* cdbEnt = cdbMgr->Get("GRP/Geometry/Data", refRun);
512  if (!cdbEnt) {
513  AliErrorF("No geometry found from %d", refRun);
514  return false;
515  }
516  // Initialize the geometry manager
517  AliInfo("Set Geometry");
518  AliGeomManager::SetGeometry(static_cast<TGeoManager*>(cdbEnt->GetObject()));
519  // Now perform mis-alignment - again based on an LHC10h run!
520  AliInfo("Misalign geometry");
521  if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",refRun,-1,-1)) {
522  AliErrorF("Failed to misalign geometry from %d", refRun);
523  return false;
524  }
525  return true;
526 }
527 //____________________________________________________________________
529 {
530  fTracklets = new TClonesArray(TrackletClassName());
531  // fTracklets->SetName(GetName());
532  TObject* obj = fTracklets;
533  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
534  AliAODHandler* ah =
535  dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
536  if (!ah) {
537  AliWarning("No AOD output handler set in analysis manager");
538  return false;
539  }
540  ah->AddBranch("TClonesArray", &obj);
541  // AliAODEvent* aod = ah->GetAOD();
542  // aod->Print();
543  // ah->GetTree()->Print();
544 
545  if (fFilterStrange > 0) {
546  fStrangeLoss = new TParameter<double>("strLoss", 0);
547  ah->AddBranch("TParameter<double>", &fStrangeLoss);
548  }
549  return true;
550 }
551 
552 //____________________________________________________________________
554 {
555  if (!InitCDB()) return false;
556  if (!InitBranch()) return false;
557 
558  // PostData(1, fContainer);
559 
560  return true;
561 }
562 //____________________________________________________________________
564 {
565  // Clear container of tracklets
566  if (!fTracklets) {
567  AliWarning("Tracklets container not initialized, init must have failed!");
568  return;
569  }
570  fTracklets->Clear();
571 
572  if (!ProcessEvent()) return;
573 
574  MarkForStore();
575 }
576 
577 //____________________________________________________________________
579 {
580  AliVEvent* event = FindEvent();
581  if (!event) return false;
582 
583 
584  TTree* clusters = FindClusters();
585  const AliVVertex* ip = FindIP(event);
586  Bool_t ret = true;
587  if (fDebug > 0) AliInfoF("Cluster tree: %p, ip: %p", clusters, ip);
588  if (!ip) return false;
589  if (clusters) {
590  if (!HasGeometry()) ret = false;
591  if (fDebug > 1) AliInfoF("After geo check: %d", ret);
592  if (ret && !HasField(event)) ret = false;
593  if (fDebug > 1) AliInfoF("After field check: %d", ret);
594  // If we have clusters, then try to reconstruct the event (again)
595  if (ret) ret = Reconstruct(clusters, ip);
596  if (fDebug > 1) AliInfoF("After reconstruction: %d", ret);
597  }
598  else
599  // Otherwise, use the stored tracklets (max Delta = 1.5)
600  ret = ProcessTracklets(true, event->GetMultiplicity());
601 
602  CleanClusters(clusters);
603  if (fDebug > 1) AliInfoF("Return value: %d", ret);
604  return ret;
605 }
606 
607 //____________________________________________________________________
609 {
610  if (!AliGeomManager::GetGeometry()) {
611  AliError("No geometry loaded, needed for reconstruction");
612  AliError("Add the AliTaskCDBconnect to the train");
613  return false;
614  }
615  return true;
616 }
617 //____________________________________________________________________
619 {
620  if (!TGeoGlobalMagField::Instance()->GetField() &&
621  !event->InitMagneticField()) {
622  AliWarning("Failed to initialize magnetic field");
623  return false;
624  }
625  return true;
626 }
627 
628 //____________________________________________________________________
630 {
631  AliVEvent* event = InputEvent();
632  if (fDebug > 0 && !event) AliWarning("No event");
633  return event;
634 }
635 
636 //____________________________________________________________________
638 {
639  AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
640  AliVEventHandler* inh = mgr->GetInputEventHandler();
641  if (!inh->IsA()->InheritsFrom(AliESDInputHandlerRP::Class())) {
642  if (fDebug > 0)
643  AliWarningF("Clusters not available via input handler of class: %s",
644  inh->ClassName());
645  return 0;
646  }
647  AliESDInputHandlerRP* rph = static_cast<AliESDInputHandlerRP*>(inh);
648  TTree* tree = rph->GetTreeR("ITS");
649  if (!tree) {
650  AliWarning("Tree of clusters (rec.points) not found");
651  return 0;
652  }
653  return FilterClusters(tree);
654 }
655 
656 //____________________________________________________________________
657 const AliVVertex* AliTrackletAODTask::FindIP(AliVEvent* event,
658  Double_t maxDispersion,
659  Double_t maxZError)
660 {
661  const AliVVertex* ip = event->GetPrimaryVertex();
662  if (ip->GetNContributors() <= 0) {
663  if (fDebug > 0)
664  AliWarning("Not enough contributors for IP");
665  return 0;
666  }
667  // If this is from the Z vertexer, do some checks
668  if (ip->IsFromVertexerZ()) {
669  // Get covariance matrix
670  Double_t covar[6];
671  ip->GetCovarianceMatrix(covar);
672  Double_t sigmaZ = TMath::Sqrt(covar[5]);
673  if (sigmaZ >= maxZError) {
674  if (fDebug)
675  AliWarningF("IPz resolution = %f >= %f", sigmaZ, maxZError);
676  return 0;
677  }
678 
679  // If this IP doesn not derive from AliVertex, don't check dispersion.
680  if (ip->IsA()->InheritsFrom(AliVertex::Class())) {
681  const AliVertex* ipv = static_cast<const AliVertex*>(ip);
682  // Dispersion is the parameter used by the vertexer for finding the IP.
683  if (ipv->GetDispersion() >= maxDispersion) {
684  if (fDebug > 0)
685  AliWarningF("IP dispersion = %f >= %f",
686  ipv->GetDispersion(), maxDispersion);
687  return 0;
688  }
689  }
690  }
691 
692 #if 0
693  // If we get here, we either have a full 3D vertex or track
694  // vertex, and we should check if it is in range
695  if (ip->GetZ() < fIPzAxis.GetXmin() || ip->GetZ() > fIPzAxis.GetXmax()) {
696  AliWarningF("IPz = %fcm out of range [%f,%f]cm",
697  ip->GetZ(), fIPzAxis.GetXmin(), fIPzAxis.GetXmax());
698  return 0;
699  }
700 #endif
701  // Good vertex, return it
702  return ip;
703 }
704 
705 //____________________________________________________________________
707 {
708 
709  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
710  AliAODHandler* ah =
711  dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
712  if (!ah) return false;
713  ah->SetFillAOD(kTRUE);
714  if (fDebug > 0) AliInfo("Storing event");
715  return true;
716 }
717 
718 //____________________________________________________________________
719 Bool_t AliTrackletAODTask::Reconstruct(TTree* clusters, const AliVVertex* ip)
720 {
721  // Make an IP array for the reconstructor
722  Float_t ipv[3];
723  ipv[0] = ip->GetX();
724  ipv[1] = ip->GetY();
725  ipv[2] = ip->GetZ();
726  if (fDebug > 0)
727  AliInfoF("Reconstructing from cluster tree %p with (%4.2f,%4.2f,%7.4f",
728  clusters, ipv[0], ipv[1], ipv[2]);
729 
730  // Make the reconstructor
731  fReco = 0;
732  AliITSMultRecBg reco;
733  reco.SetCreateClustersCopy (true);
734  reco.SetScaleDThetaBySin2T (fScaleDTheta);
735  reco.SetNStdDev (fMaxDelta);
736  reco.SetPhiWindow (fDPhiWindow);
737  reco.SetThetaWindow (fDThetaWindow);
738  reco.SetPhiShift (fDPhiShift);
739  reco.SetRemoveClustersFromOverlaps(fPhiOverlapCut > 0 ||
740  fZEtaOverlapCut > 0);
741  reco.SetPhiOverlapCut (fPhiOverlapCut);
742  reco.SetZetaOverlapCut (fZEtaOverlapCut);
743  reco.SetHistOn (false);
744  reco.SetRecType (AliITSMultRecBg::kData);
745 
746  // Run normal reconstruction
747  reco.Run(clusters, ipv);
748 
749  // Save pointer to reco object - for MC combinatorics
750  fReco = &reco;
751  // And fill results into output branch
752  if (!ProcessTracklets(true, reco.GetMultiplicity())) {
753  AliWarning("Process tracklets (normal) failed");
754  return false;
755  }
756  fReco = 0;
757 
758  // Run again, but in injection mode
759  reco.SetRecType(AliITSMultRecBg::kBgInj);
760  reco.Run(clusters, ipv);
761 
762  // And fill results into output branch
763  if(!ProcessTracklets(false, reco.GetMultiplicity())) {
764  AliWarning("Process tracklets (injection) failed");
765  return false;
766  }
767  if (fDebug > 1) AliInfo("Returning true");
768  return true;
769 }
770 
771 //____________________________________________________________________
773  AliVMultiplicity* vmult)
774 {
775  // Now loop over all tracklets. Since this method is only called
776  // once per "reconstruction" pass, we save a lot of computing time
777  // since we loop over the tracklets of each reconstruction exactly
778  // once.
779  // if (!mult->IsA()->InheritsFrom(AliMultiplicity::Class())) return false;
780  AliMultiplicity* mult = static_cast<AliMultiplicity*>(vmult);
781  Int_t nTracklets = mult->GetNumberOfTracklets();
782  for (Int_t trackletNo = 0; trackletNo < nTracklets; trackletNo++) {
783  if (!ProcessTracklet(normal,mult,trackletNo)) return false;
784  }
785  return true;
786 }
787 
788 //____________________________________________________________________
790 {
791  if (!fTracklets) return 0;
792  Int_t n = fTracklets->GetEntries();
793  return new((*fTracklets)[n]) AliAODTracklet;
794 }
795 //____________________________________________________________________
798  AliMultiplicity* mult,
799  Int_t no)
800 {
801  Double_t theta = mult->GetTheta (no);
802  Double_t phi = mult->GetPhi (no);
803  Double_t dTheta = mult->GetDeltaTheta(no);
804  Double_t dPhi = mult->GetDeltaPhi (no);
805  Double_t delta = mult->CalcDist (no);
806  AliAODTracklet* tracklet = MakeTracklet(normal);
807  tracklet->SetTheta (theta);
808  tracklet->SetPhi (phi);
809  tracklet->SetDTheta(dTheta);
810  tracklet->SetDPhi (dPhi);
811  tracklet->SetDelta (delta);
812  if (!normal) tracklet->SetInjection();
813 
814  return tracklet;
815 }
816 
817 /*********************************************************************
818  *
819  * Code for processing simulations
820  */
821 #ifndef __CINT__
822 #include <AliStack.h>
823 #include <AliMCEvent.h>
824 #include <AliGenEventHeader.h>
825 #include <TParticle.h>
826 #include <TParticlePDG.h>
827 #include <TDatabasePDG.h>
828 #include <AliITSgeomTGeo.h>
829 #include <TRandom.h>
830 #include <TROOT.h>
831 #else
832 class TParticlePDG;
833 class TParticle;
834 class AliStack;
835 #endif
836 //====================================================================
843 {
844 public:
854  AliTrackletAODMCTask(const char* name)
855  : AliTrackletAODTask(name)
856  {}
863  : AliTrackletAODTask(other)
864  {}
873  {
874  return *this;
875  }
876 protected:
881  /*
882  * Initialize the worker
883  */
884  Bool_t WorkerInit();
890  virtual TTree* FilterClusters(TTree* t);
901  virtual void FilterClustersRandom(TTree* t, TTree* copy);
913  virtual void FilterClustersTrack(TTree* t, TTree* copy);
919  virtual void CleanClusters(TTree*& t);
926  virtual Bool_t ProcessEvent();
933  /* @} */
943  virtual const char* GetCDBReferenceURL() const
944  {
945  return "alien://Folder=/alice/simulation/2008/v4-15-Release/Residual";
946  }
947  /* @} */
957  virtual const char* TrackletClassName() const { return "AliAODMCTracklet"; }
965  virtual AliAODTracklet* MakeTracklet(Bool_t normal);
975  virtual AliAODTracklet* ProcessTracklet(Bool_t normal,
976  AliMultiplicity* mult,
977  Int_t no);
978  /* @} */
992  Int_t FindPrimaryParentID(Int_t label) const;
1001  TParticle* FindPrimaryParent(Int_t label) const;
1010  Int_t FindParent(Int_t label) const;
1025  Int_t FindParents(Int_t label, TArrayI& fill, Int_t offset) const;
1036  Int_t CommonParent(Int_t label, const TArrayI& fill) const;
1037  /* @} */
1038 
1039 
1041 };
1042 
1043 #define ADDP(P,M,L,Q) \
1044  db->AddParticle("p" #P, "p" #P,M,false,L,Q,"Baryon",P); \
1045  db->AddAntiParticle("p" #P "_bar", -P)
1046 
1047 //____________________________________________________________________
1049 {
1050  if (!AliTrackletAODTask::WorkerInit()) return false;
1051 
1052  // Register additional particles from EPOS-LHC
1053  TDatabasePDG* db = TDatabasePDG::Instance();
1054  db->GetParticle(2212);
1055  db->AddParticle("deut","deut",1.876560,false,0.000000,3,"Baryon",1000010020);
1056  db->AddAntiParticle("deut_bar",-1000010020);
1057  db->AddParticle("trit","trit",2.816700,false,0.000000,3,"Baryon",1000010030);
1058  db->AddAntiParticle("trit_bar",-1000010030);
1059  db->AddParticle("alph","alph",3.755000,false,0.000000,6,"Baryon",1000020040);
1060  db->AddAntiParticle("alph_bar",-1000020040);
1061  ADDP(32224, 1.524000,0.145000,6);
1062  ADDP(12224, 1.816000,0.350000,6);
1063  ADDP(12222, 2.108000,0.350000,6);
1064  ADDP(12212, 1.525720,0.145000,3);
1065  ADDP(2124, 1.819440,0.150000,3);
1066  ADDP(32214, 2.113160,0.150000,3);
1067  ADDP(32212, 2.406880,0.145000,3);
1068  ADDP(12214, 2.700600,0.300000,3);
1069  ADDP(22124, 2.994320,0.150000,3);
1070  ADDP(12122, 3.288040,0.300000,3);
1071  ADDP(13222, 1.575200,0.080000,3);
1072  ADDP(23222, 1.768100,0.100000,3);
1073  ADDP(13226, 1.961000,0.170000,3);
1074  ADDP(12112, 1.524430,0.350000,0);
1075  ADDP(1214, 1.816860,0.150000,0);
1076  ADDP(32114, 2.109290,0.150000,0);
1077  ADDP(32112, 2.401720,0.145000,0);
1078  ADDP(12114, 2.694150,0.300000,0);
1079  ADDP(21214, 2.986579,0.150000,0);
1080  ADDP(11212, 3.279010,0.300000,0);
1081  ADDP(13122, 1.761000,0.040000,0);
1082  ADDP(3124, 1.950500,0.016000,0);
1083  ADDP(23122, 2.140000,0.090000,0);
1084  ADDP(13212, 2.329500,0.080000,0);
1085  ADDP(23212, 2.519000,0.100000,0);
1086  ADDP(43122, 2.708500,0.145000,0);
1087  ADDP(13216, 2.898000,0.170000,0);
1088  ADDP(31114, 1.524000,0.145000,-3);
1089  ADDP(11114, 1.816000,0.350000,-3);
1090  ADDP(11112, 2.108000,0.350000,-3);
1091  ADDP(13112, 1.577600,0.080000,-3);
1092  ADDP(23112, 1.767700,0.100000,-3);
1093  ADDP(13116, 1.957800,0.170000,-3);
1094  ADDP(9900110,0.000000,0.000000,0);
1095  ADDP(9900210,0.000000,0.000000,0);
1096  ADDP(9900220,0.000000,0.000000,0);
1097  ADDP(9900330,0.000000,0.000000,0);
1098  ADDP(9900440,0.000000,0.000000,0);
1099  ADDP(9902210,0.000000,0.000000,0);
1100  ADDP(9902110,0.000000,0.000000,0);
1101  ADDP(88, 0.000000,0.000000,0);
1102  ADDP(90, 0.000000,0.000000,0);
1103  ADDP(990, 0.000000,0.000000,0);
1104  ADDP(99999, 0.000000,0.000000,0);
1105 
1106  return true;
1107 }
1108 
1109 namespace {
1110  const Double_t k0s = 1.52233299626516083e+00; // 310 - K^0_S weight
1111  const Double_t kpm = (1.43744204476109627e+00*
1112  9.82150320171071400e-01); // 321 - K^{+/-}
1113  const Double_t lam = 2.75002089647900005e+00; // 3122 - lambda
1114  const Double_t sig = 2.75002089647899961e+00; // 3212 - sigma
1115  const Double_t xi = 3.24109605656453548e+00; // 3322 - Xi
1116 
1118  {
1119  switch (TMath::Abs(pdg)) {
1120  case 310: return k0s;
1121  case 321: return kpm;
1122  case 3122: return lam;
1123  case 3212: return sig;
1124  case 3322: return xi;
1125  }
1126  return 1;
1127  }
1128 
1129  Bool_t KeepIt(Double_t weight)
1130  {
1131  if (weight <= 1) return true;
1132  Double_t chance = 1 - 1 / weight; // chance is 1 minus inverse
1133  return (gRandom->Uniform() >= chance);
1134  }
1135 }
1136 //____________________________________________________________________
1138 {
1139  if (fStrangeLoss) fStrangeLoss->SetVal(0);
1140  if (!t || fFilterStrange <= 0) {
1141  if (fDebug > 0) AliInfo("Returning original cluster tree");
1142  return t;
1143  }
1144 
1145  TDirectory* savDir = gDirectory;
1146  gDirectory = gROOT;
1147  TTree* copy = new TTree("TreeR", "TreeR");
1148  copy->SetAutoFlush(0);// Keep in memory
1149 
1150  if (fFilterStrange == 1) FilterClustersRandom(t, copy);
1151  else FilterClustersTrack(t, copy);
1152 
1153  savDir->cd();
1154  if (fDebug > 0)
1155  AliInfoF("Returning reduced cluster tree %p (original %p)", copy, t);
1156  return copy;
1157 }
1158 
1159 
1160 //____________________________________________________________________
1162 {
1163  TClonesArray* in = new TClonesArray("AliITSRecPoint");
1164  TClonesArray* out = new TClonesArray("AliITSRecPoint");
1165  Int_t outN = 0;
1166  copy->Branch("ITSRecPoints", &out);
1167  t->SetBranchAddress("ITSRecPoints", &in);
1168 
1169  // Printf("Filtering clusters from strange particles");
1170 
1171  Int_t min1 = AliITSgeomTGeo::GetModuleIndex(1,1,1);
1172  Int_t max1 = AliITSgeomTGeo::GetModuleIndex(2,1,1);
1173  Int_t min2 = AliITSgeomTGeo::GetModuleIndex(2,1,1);
1174  Int_t max2 = AliITSgeomTGeo::GetModuleIndex(3,1,1);
1175  Int_t nTotal = 0;
1176  Int_t nKept = 0;
1177 
1178  // Loop over the modules of the SPD
1179  for (Int_t i = 0; i < max2; i++) {
1180  in->Clear();
1181  out->Clear();
1182  outN = 0;
1183 
1184  // Read in module data
1185  t->GetEntry(i);
1186 
1187  // Loop over all clusters in the module
1188  Int_t inN = in->GetEntries();
1189  for (Int_t j = 0; j < inN; j++) {
1190  AliITSRecPoint* inCl = static_cast<AliITSRecPoint*>(in->At(j));
1191  if (!inCl) continue;
1192 
1193  // Loop over labels in the module
1194  Double_t weight = 1;
1195  for (Int_t k = 0; k < 3; k++) {
1196  Int_t label = inCl->GetLabel(k);
1197  if (label <= 0) continue;
1198 
1199  // Check primary parent particle type
1200  TParticle* parent = FindPrimaryParent(label);
1201  if (!parent) continue;
1202 
1203  weight *= GetWeight(parent->GetPdgCode());
1204  if (weight > 1)
1205  Printf("Cluster from %+5d: %7.5f", parent->GetPdgCode(), weight);
1206  // Printf("Parent %d is a K^0_S: %d", k+1, parent->GetPdgCode());
1207  // Randomly remove clusters from the right kind of primary
1208  // parent particles
1209  }
1210  if (weight > 1) {
1211  nTotal++;
1212  if (!KeepIt(weight)) continue;
1213  nKept++;
1214  }
1215  new ((*out)[outN++]) AliITSRecPoint(*inCl);
1216  }
1217  // Printf("Kept %d out of %d clusters", outN, inN);
1218  copy->Fill();
1219  }
1220  if (nTotal > 0) {
1221  Double_t loss = 1-Float_t(nKept)/nTotal;
1222  Printf("Kept %d out of %d clusters from strange primaries (%4.1f%%)",
1223  nKept, nTotal, 100.*loss);
1224  if (fStrangeLoss) fStrangeLoss->SetVal(loss);
1225  }
1226  delete in;
1227  delete out;
1228 }
1229 
1230 //____________________________________________________________________
1232 {
1233  TClonesArray* in = new TClonesArray("AliITSRecPoint");
1234  TClonesArray* out = new TClonesArray("AliITSRecPoint");
1235  Int_t outN = 0;
1236  copy->Branch("ITSRecPoints", &out);
1237  t->SetBranchAddress("ITSRecPoints", &in);
1238 
1239  // Loop over the stack
1240  Int_t nTotal = 0;
1241  Int_t nKept = 0;
1242  AliStack* stack = MCEvent()->Stack();
1243  Int_t nTracks = stack->GetNtrack();
1244  TBits kept(nTracks);
1245  TBits seen(nTracks);
1246  kept.ResetAllBits(true);
1247  seen.ResetAllBits(false);
1248  for (Int_t trackNo = nTracks; trackNo--; ) {
1249  // Get the primary parent identifier
1250  Int_t parent = FindPrimaryParentID(trackNo);
1251  if (parent < 0) continue;
1252  // Check if we've seen this parent already
1253  if (seen.TestBitNumber(parent)) continue;
1254  seen.SetBitNumber(parent,true);
1255 
1256  // Get the primary parent track and weight
1257  TParticle* par = stack->Particle(parent);
1258  Double_t weight = GetWeight(par->GetPdgCode());
1259  Bool_t keep = true;
1260  if (weight > 1) {
1261  nTotal++;
1262  keep = KeepIt(weight);
1263  if (keep) nKept++;
1264  }
1265  if (fDebug > 3)
1266  Printf("Primary parent %6d from a %6d %6s (%f)",
1267  parent, par->GetPdgCode(), keep ? "kept" : "marked", weight);
1268  kept.SetBitNumber(parent, keep);
1269  par->SetWeight(keep ? 1 : weight);
1270  }
1271  // At this point, we have investigated all relevant primary
1272  // particles and thrown a dice to see if we should remove clusters
1273  // that correspond to these primary particles.
1274  if (fDebug > 0 && nTotal > 0) {
1275  Double_t loss = 1-Float_t(nKept)/nTotal;
1276  Printf("Kept %d out of %d strange primaries (%4.1f%%)",
1277  nKept, nTotal, 100.*loss);
1278  if (fStrangeLoss) fStrangeLoss->SetVal(loss);
1279  }
1280 
1281  // Next thing is to actually remove the clusters
1282  // Printf("Filtering clusters from strange particles");
1283  Int_t min1 = AliITSgeomTGeo::GetModuleIndex(1,1,1);
1284  Int_t max1 = AliITSgeomTGeo::GetModuleIndex(2,1,1);
1285  Int_t min2 = AliITSgeomTGeo::GetModuleIndex(2,1,1);
1286  Int_t max2 = AliITSgeomTGeo::GetModuleIndex(3,1,1);
1287  Int_t inT = 0;
1288  Int_t outT = 0;
1289  // Loop over the modules of the SPD
1290  for (Int_t i = 0; i < max2; i++) {
1291  in->Clear();
1292  out->Clear();
1293  outN = 0;
1294 
1295  // Read in module data
1296  t->GetEntry(i);
1297 
1298  // Loop over all clusters in the module
1299  Int_t inN = in->GetEntries();
1300  inT += inN;
1301  for (Int_t j = 0; j < inN; j++) {
1302  AliITSRecPoint* inCl = static_cast<AliITSRecPoint*>(in->At(j));
1303  if (!inCl) continue;
1304 
1305  // Loop over labels in the module
1306  Bool_t toRemove = false;
1307  for (Int_t k = 0; k < 3; k++) {
1308  Int_t label = inCl->GetLabel(k);
1309  if (label <= 0) continue;
1310 
1311  // Check primary parent particle type
1312  Int_t parent = FindPrimaryParentID(label);
1313  if (parent < 0) continue;
1314 
1315  if (!kept.TestBitNumber(parent)) toRemove = true;
1316  if (fDebug > 3) {
1317  Printf("Cluster %6d from parent %6d (%7.5f) of type %6d %s",
1318  j, parent,
1319  stack->Particle(parent)->GetPdgCode(),
1320  stack->Particle(parent)->GetWeight(),
1321  toRemove ? "removed" : "kept");
1322  }
1323  }
1324  // We've now looked at all primaries for this cluster, and if
1325  // any of these primaries was flagged for removal, then we
1326  // should remove the cluster.
1327  if (toRemove) continue;
1328  new ((*out)[outN++]) AliITSRecPoint(*inCl);
1329  }
1330  // Printf("Kept %d out of %d clusters", outN, inN);
1331  outT += outN;
1332  copy->Fill();
1333  }
1334  if (fDebug > 0)
1335  Printf("Wrote out %6d out of %6d clusters (%4.1f%%)",
1336  outT, inT, (inT > 0 ? 100*float(inT-outT)/inT : 100));
1337  delete in;
1338  delete out;
1339 }
1340 
1341 //____________________________________________________________________
1343 {
1344  if (!t || fFilterStrange <= 0) return;
1345 
1346  delete t;
1347  t = 0;
1348 }
1349 
1350 //____________________________________________________________________
1352 {
1353  if (!AliTrackletAODTask::ProcessEvent()) return false;
1354 
1355  // create generated "tracklets"
1356  if (!ProcessGenerated()) return false;
1357 
1358  return true;
1359 }
1360 
1361 //____________________________________________________________________
1363 {
1364  if (fDebug > 0) AliInfo("Processing generated particles");
1365  AliStack* stack = MCEvent()->Stack();
1366  Int_t nTracks = stack->GetNtrack();
1367  for (Int_t trackNo = nTracks; trackNo--; ) {
1368  if (!stack->IsPhysicalPrimary(trackNo)) {
1369  // AliWarningF("Track # %6d not a primary", trackNo);
1370  // Not a primary, go on
1371  continue;
1372  }
1373  // Get the particle
1374  TParticle* particle = stack->Particle(trackNo);
1375  if (!particle) {
1376  AliWarningF("No particle found for track # %d", trackNo);
1377  continue;
1378  }
1379  TParticlePDG* pdg = particle->GetPDG();
1380  if (!pdg) {
1381  AliWarningF("Unknown PDG code: %d", particle->GetPdgCode());
1382  // continue; // Unknown particle
1383  }
1384  // if (pdg->Charge() == 0) {
1385  // Uncharged - don't care
1386  // continue;
1387  // }
1388 
1389  // Get theta
1390  Double_t theta = particle->Theta();
1391  // Check for beam-like particle
1392  if (theta < 1e-6 || TMath::Abs(theta-TMath::Pi()) < 1e-6) {
1393  if (fDebug > 0)
1394  AliWarningF("Track # %6d is beam-like (%f)", trackNo,
1395  TMath::RadToDeg()*theta);
1396  continue;
1397  }
1398  Double_t eta = -TMath::Log(TMath::ATan(theta/2));
1399  // If the pseudorapidity is way beyond the SPD acceptance, do not
1400  // write down this generated "tracklet" - to save space.
1401  if (TMath::Abs(eta) > 3) continue;
1402 
1403  Double_t phi = particle->Phi();
1404  AliAODTracklet* tracklet = MakeTracklet(false);
1405  AliAODMCTracklet* mc = static_cast<AliAODMCTracklet*>(tracklet);
1406  mc->SetGenerated();
1407  mc->SetTheta(theta);
1408  mc->SetPhi(phi);
1409  mc->SetDTheta(0);
1410  mc->SetDPhi (0);
1411  mc->SetDelta (0);
1412  mc->SetParentPdg(particle->GetPdgCode());
1413  mc->SetParentPt (particle->Pt());
1414  if (particle->GetWeight() > 1) mc->SetSuppressed();
1415  if (pdg && pdg->Charge() == 0) mc->SetNeutral();
1416  }
1417  if (fDebug > 2) AliInfo("Returning true from generated");
1418  return true;
1419 }
1420 
1421 //____________________________________________________________________
1423 {
1424  if (!fTracklets) return 0;
1425  Int_t n = fTracklets->GetEntries();
1426  return new((*fTracklets)[n]) AliAODMCTracklet;
1427 }
1428 //____________________________________________________________________
1431  AliMultiplicity* mult,
1432  Int_t no)
1433 {
1434  AliAODTracklet* tracklet =
1435  AliTrackletAODTask::ProcessTracklet(normal,mult,no);
1436  if (!normal) return tracklet;
1437 
1438  AliAODMCTracklet* mc = static_cast<AliAODMCTracklet*>(tracklet);
1439  Int_t label0 = mult->GetLabel(no, 0);
1440  Int_t label1 = mult->GetLabel(no, 1);
1441  TParticle* parent0 = FindPrimaryParent(label0);
1442  if (parent0) {
1443  mc->SetParentPdg(parent0->GetPdgCode());
1444  mc->SetParentPt (parent0->Pt());
1445  }
1446  if (label0 != label1) {
1447  TParticle* parent1 = FindPrimaryParent(label1);
1448  if (parent1) {
1449  mc->SetParentPdg(parent1->GetPdgCode(), true);
1450  mc->SetParentPt (parent1->Pt(), true);
1451  }
1452  mc->SetCombinatorics();
1453  // Here, we could track back in the cluster labels to see if we
1454  // have the same ultimate mother of both clusters. Since this
1455  // isn't really used, we do not do that
1456  if (fReco) {
1457  TArrayI parents(50); // At most 50 levels deep
1458  // Tracklet parameters
1459  Float_t* ftrack = fReco->GetTracklet(no);
1460  // Cluster identifiers
1461  Int_t clus0 = Int_t(ftrack[AliITSMultReconstructor::kClID1]);
1462  Int_t clus1 = Int_t(ftrack[AliITSMultReconstructor::kClID2]);
1463  // Cluster labelsx
1464  Float_t* fclus0 = (fReco->GetClusterOfLayer(0,clus0) +
1465  AliITSMultReconstructor::kClMC0);
1466  Float_t* fclus1 = (fReco->GetClusterOfLayer(1,clus1) +
1467  AliITSMultReconstructor::kClMC0);
1468  // Loop over three inner layer cluster labels
1469  Int_t offset = 0;
1470  for (Int_t i = 0; i < 3; i++)
1471  offset = FindParents(Int_t(fclus0[i]), parents, offset);
1472  // Loop over three outer layer cluster labels
1473  Bool_t distinct = true;
1474  for (Int_t i = 0; i < 3; i++) {
1475  Float_t flbl = fclus1[i];
1476  // Be careful not to get overflows on conversions
1477  if (flbl > nextafter(INT_MAX, 0) || flbl < nextafter(INT_MIN, 0))
1478  continue; // Would overflow
1479  if (CommonParent(Int_t(flbl), parents)) {
1480  distinct = false;
1481  // We break out as soon as we find a common parent.
1482  break;
1483  }
1484  } // loop over outer layer
1485  if (distinct) mc->SetDistinct();
1486  } // if (fReco)
1487  }
1488  else {
1489  if (!MCEvent()->Stack()->IsPhysicalPrimary(label0))
1490  mc->SetSecondary();
1491  }
1492  return tracklet;
1493 }
1494 
1495 //____________________________________________________________________
1497 {
1498  AliStack* stack = MCEvent()->Stack();
1499  Int_t nTracks = stack->GetNtrack();
1500  TParticle* particle = stack->Particle(label);
1501  if (!particle) return -1;
1502  Int_t ret = particle->GetFirstMother();
1503  if (ret > nTracks) return -1;
1504  return ret;
1505 }
1506 
1507 //____________________________________________________________________
1509  TArrayI& fill,
1510  Int_t offset) const
1511 {
1512  // If offset is negative, then we don't really store the labels, but
1513  // only check against those already stored.
1514  Int_t i = offset;
1515  Int_t lbl = label;
1516  while (i < fill.GetSize()-1 && lbl >= 1) {
1517  fill[i] = lbl;
1518  i++;
1519  lbl = FindParent(lbl);
1520  }
1521  // If we get here, and we're checking, that means we did not find a
1522  // common ancestor, so we return 0 (false). If we're not checking,
1523  // then we return the number of elements assigned in the passed
1524  // cache.
1525  return i;
1526 }
1527 
1528 //____________________________________________________________________
1530  const TArrayI& fill) const
1531 {
1532  Int_t i = 0;
1533  Int_t lbl = label;
1534  while (i < fill.GetSize()-1 && lbl >= 1) {
1535  // If we're checking, just see if we have the label in the list
1536  // already. If we do, then return the index
1537  for (Int_t j = 0; j < fill.GetSize(); j++) {
1538  if (fill[j] == lbl) return j;
1539  }
1540  i++;
1541  lbl = FindParent(lbl);
1542  }
1543  // If we get here, and we're checking, that means we did not find a
1544  // common ancestor, so we return 0 (false). If we're not checking,
1545  // then we return the number of elements assigned in the passed
1546  // cache.
1547  return false;
1548 }
1549 
1550 
1551 //____________________________________________________________________
1553 {
1554  Int_t parent = FindPrimaryParentID(label);
1555  if (parent < 0) return 0;
1556  return MCEvent()->Stack()->Particle(parent);
1557 }
1558 //____________________________________________________________________
1560 {
1561  AliStack* stack = MCEvent()->Stack();
1562  Int_t nTracks = stack->GetNtrack();
1563  Int_t trackNo = label;
1564  if (trackNo > nTracks || trackNo <= 0) return 0;
1565  TParticle* particle = stack->Particle(label);
1566  while (!stack->IsPhysicalPrimary(trackNo)) {
1567  trackNo = particle->GetFirstMother();
1568  // If we have hit the top
1569  if (trackNo < 0) return 0;
1570  // Partice first next iteration
1571  particle = stack->Particle(trackNo);
1572  }
1573  return trackNo;
1574 }
1575 //====================================================================
1577 {
1578  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1579  if (!mgr) {
1580  ::Error("Create","No analysis manager to connect to.");
1581  return 0;
1582  }
1583  AliAODHandler* ah =
1584  dynamic_cast<AliAODHandler*>(mgr->GetOutputEventHandler());
1585  if (!ah) {
1586  ::Error("Create","No AOD output handler!");
1587  return 0;
1588  }
1589  Bool_t mc = mgr->GetMCtruthEventHandler() != 0;
1590  AliTrackletAODTask* ret = 0;
1591  if (mc) ret = new AliTrackletAODMCTask("MidRapidityMC");
1592  else ret = new AliTrackletAODTask("MidRapidity");
1593  if (ret) ret->Connect();
1594 
1595  return ret;
1596 }
1597 
1598 
1599 //====================================================================
1600 
1601 //
1602 // EOF
1603 //
1604 
1605 
Int_t pdg
void SetPhiOverlapCut(Double_t x=0.005)
void SetPhi(Real_t x)
double Double_t
Definition: External.C:58
void UserExec(Option_t *)
Double_t GetWeight(Int_t averageoption, Double_t pt, TH1D *hRaa, Double_t raaSystLow, Double_t raaSystHigh, Double_t ppSystRawYieldCutVar, Double_t ppSystRawYieldCutVarPid, Double_t ABSystRawYieldCutVar, Double_t ABSystRawYieldCutVarPid)
Real_t GetTheta() const
virtual const char * GetCDBReferenceURL() const
Int_t FindPrimaryParentID(Int_t label) const
virtual AliAODTracklet * MakeTracklet(Bool_t normal)
virtual Bool_t WorkerInit()
void SetCombinatorics()
Int_t FindParent(Int_t label) const
void SetParentPdg(Short_t pdg, Bool_t second=false)
AliTrackletAODTask & operator=(const AliTrackletAODTask &other)
void SetTheta(Real_t x)
AliTrackletAODMCTask(const char *name)
AliITSMultRecBg * fReco
void SetDPhiWindow(Double_t x=0.06)
void SetDTheta(Real_t x)
TRandom * gRandom
void SetDPhiShift(Double_t x=0.0045)
virtual void CleanClusters(TTree *&t)
ClassDef(AliTrackletAODTask, 1)
virtual Bool_t ProcessEvent()
TString kData
Declare data MC or deltaAOD.
Tracklet AOD object.
virtual TTree * FilterClusters(TTree *t)
TParticle * FindPrimaryParent(Int_t label) const
virtual const char * TrackletClassName() const
static AliTrackletAODTask * Create()
int Int_t
Definition: External.C:63
Int_t CommonParent(Int_t label, const TArrayI &fill) const
virtual AliAODTracklet * ProcessTracklet(Bool_t normal, AliMultiplicity *mult, Int_t no)
float Float_t
Definition: External.C:68
void SetZEtaOverlapCut(Double_t x=0.05)
Bool_t ProcessTracklets(Bool_t normal, AliVMultiplicity *mult)
TParameter< double > * fStrangeLoss
virtual Bool_t ProcessEvent()
Int_t mode
Definition: anaM.C:40
virtual Int_t GetCDBReferenceRun() const
void SetParentPt(Real_t pt, Bool_t second=false)
virtual AliAODTracklet * MakeTracklet(Bool_t normal)
void SetScaleDTheta(Bool_t x=false)
void SetDThetaWindow(Double_t x=0.025)
const AliVVertex * FindIP(AliVEvent *event, Double_t maxDispersion=0.04, Double_t maxZError=0.25)
virtual void CleanClusters(TTree *&t)
Bool_t HasField(AliVEvent *event)
TClonesArray * fTracklets
void SetMaxDelta(Double_t x=25)
virtual AliAODTracklet * ProcessTracklet(Bool_t normal, AliMultiplicity *mult, Int_t no)
virtual TTree * FilterClusters(TTree *t)
ClassDef(AliTrackletAODMCTask, 1)
void Print(Option_t *) const
virtual void FilterClustersRandom(TTree *t, TTree *copy)
virtual const char * GetCDBReferenceURL() const
AliTrackletAODMCTask(const AliTrackletAODMCTask &other)
virtual void SetFilterStrange(Int_t mode)
const char Option_t
Definition: External.C:48
Bool_t Reconstruct(TTree *clusters, const AliVVertex *ip)
AliTrackletAODMCTask & operator=(const AliTrackletAODMCTask &other)
bool Bool_t
Definition: External.C:53
void SetDelta(Real_t x)
void SetDPhi(Real_t x)
virtual void FilterClustersTrack(TTree *t, TTree *copy)
Int_t FindParents(Int_t label, TArrayI &fill, Int_t offset) const
#define ADDP(P, M, L, Q)
virtual const char * TrackletClassName() const
void Terminate(Option_t *)