AliPhysics  d497547 (d497547)
 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 #include "AliTrackletAODUtils.C"
14 #ifndef __CINT__
15 #include "AliAODTracklet.C"
16 #include "AliTrackletWeights.C"
17 #include "AliVVertex.h"
18 #include <AliVEvent.h>
19 #include <AliAODEvent.h>
20 #include <AliVVertex.h>
21 #include <AliVertex.h>
22 #include <AliVMultiplicity.h>
23 #include <AliMultiplicity.h>
24 #include <AliAnalysisManager.h>
25 #include <AliVEventHandler.h>
26 #include <AliESDInputHandlerRP.h>
27 #include <AliITSMultRecBg.h>
28 #include <AliCDBManager.h>
29 #include <AliCDBEntry.h>
30 #include <AliCDBPath.h>
31 #include <AliCDBId.h>
32 #include <AliGeomManager.h>
33 #include <AliAODHandler.h>
34 #include <TGeoManager.h>
35 #include <TGeoGlobalMagField.h>
36 #include <TH1.h>
37 #include <TProfile2D.h>
38 #include <TUrl.h>
39 #include <climits>
40 #else
43 class AliAODTracklet;
44 class AliVEvent;
45 class AliVVertex;
46 class AliMultiplicity;
47 class AliVMultiplicity;
48 class AliAnalysisDataContainer;
49 class AliITSMultRecBg;
50 class AliGeomManager; // Auto-load
51 class AliCDBManager; // Auto-load
52 class AliMultSelection; // Auto-load
53 class TClonesArray;
54 class TTree;
55 class TGeoGlobalMagField; // Auto-load
56 class TGeoManager; // Auto-load
57 class TParticle;
58 class TH1;
59 class TProfile2D;
60 #endif
61 
62 
63 //====================================================================
70 {
71 public:
75  enum EStatus {
83  kIP,
94  };
104  AliTrackletAODTask(const char* name);
122  void Print(Option_t*) const;
128  static AliTrackletAODTask* Create(const char* weights="");
141  void UserExec(Option_t*);
146  void FinishTaskOutput() { /*WorkerFinalize();*/ }
151  void Terminate(Option_t*) {}
157  Bool_t Connect();
158  /* @} */
168  void SetScaleDTheta(Bool_t x=false) { fScaleDTheta = x; }
174  void SetDPhiShift(Double_t x=0.0045) { fDPhiShift = x; }
180  void SetMaxDelta(Double_t x=25) { fMaxDelta = x; }
186  void SetDThetaWindow(Double_t x=0.025) { fDThetaWindow = x; }
192  void SetDPhiWindow(Double_t x=0.06) { fDPhiWindow = x; }
198  void SetPhiOverlapCut(Double_t x=0.005) { fPhiOverlapCut = x; }
214  virtual void SetFilterMode(Int_t mode) { fFilterMode = mode; }
221  /* @} */
222 protected:
230  virtual Bool_t WorkerInit();
236  Bool_t InitCDB();
242  virtual Int_t GetCDBReferenceRun() const { return 137161; }
248  virtual const char* GetCDBReferenceURL() const
249  {
250  return "alien://Folder=/alice/data/2010/OCDB";
251  }
257  Bool_t InitBranch();
268  virtual Bool_t ProcessEvent();
280  Bool_t HasField(AliVEvent* event);
286  TTree* FindClusters();
292  virtual TTree* FilterClusters(TTree* t) { return t; }
298  virtual void CleanClusters(TTree*& t) {};
304  AliVEvent* FindEvent();
314  const AliVVertex* FindIP(AliVEvent* event,
315  Double_t maxDispersion=0.04,
316  Double_t maxZError=0.25);
323  Bool_t Reconstruct(TTree* clusters, const AliVVertex* ip);
330  /* @} */
340  virtual const char* TrackletClassName() const { return "AliAODTracklet"; }
348  virtual AliAODTracklet* MakeTracklet(Bool_t normal);
357  Bool_t ProcessTracklets(Bool_t normal, AliVMultiplicity* mult);
367  virtual AliAODTracklet* ProcessTracklet(Bool_t normal,
368  AliMultiplicity* mult,
369  Int_t no);
370  /* @} */
371 
373  TClonesArray* fTracklets;
374  /* Output container */
378 
380 
395  AliITSMultRecBg* fReco;
396 
400 
401  ClassDef(AliTrackletAODTask,1);
402 };
403 //____________________________________________________________________
405  : AliAnalysisTaskSE(),
406  fTracklets(0),
407  fContainer(0),
408  fStatus(0),
409  fNClustersVsNTracklets(0),
410  fScaleDTheta(true),
411  fMaxDelta(25),
412  fDPhiShift(0.0045),
413  fDThetaWindow(0.025),
414  fDPhiWindow(0.06),
415  fPhiOverlapCut(0.005),
416  fZEtaOverlapCut(0.05),
417  fReco(0),
418  fFilterMode(0),
419  fStrangeLoss(0)
420 {}
421 //____________________________________________________________________
423  : AliAnalysisTaskSE("tracklets"),
424  fTracklets(0),
425  fContainer(0),
426  fStatus(0),
427  fNClustersVsNTracklets(0),
428  fScaleDTheta(true),
429  fMaxDelta(25),
430  fDPhiShift(0.0045),
431  fDThetaWindow(0.025),
432  fDPhiWindow(0.06),
433  fPhiOverlapCut(0.005),
434  fZEtaOverlapCut(0.05),
435  fReco(0),
436  fFilterMode(0),
437  fStrangeLoss(0)
438 {
439  DefineOutput(1,TList::Class());
440 }
441 //____________________________________________________________________
443  : AliAnalysisTaskSE(other),
444  fTracklets(0),
445  fContainer(0),
446  fStatus(0),
447  fNClustersVsNTracklets(0),
448  fScaleDTheta(other.fScaleDTheta),
449  fMaxDelta(other.fMaxDelta),
450  fDPhiShift(other.fDPhiShift),
451  fDThetaWindow(other.fDThetaWindow),
452  fDPhiWindow(other.fDPhiWindow),
453  fPhiOverlapCut(other.fPhiOverlapCut),
454  fZEtaOverlapCut(other.fZEtaOverlapCut),
455  fReco(0),
456  fFilterMode(other.fFilterMode),
457  fStrangeLoss(0)
458 {}
459 //____________________________________________________________________
461 {
462  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
463  if (!mgr) {
464  AliError("No analysis manager to connect to.");
465  return false;
466  }
467  AliAODHandler* ah =
468  dynamic_cast<AliAODHandler*>(mgr->GetOutputEventHandler());
469  if (!ah) {
470  AliError("No AOD output handler!");
471  return false;
472  }
473 
474  // Add to the manager
475  mgr->AddTask(this);
476 
477  AliAnalysisDataContainer* sumCon =
478  mgr->CreateContainer(Form("%sSums", GetName()), TList::Class(),
479  AliAnalysisManager::kOutputContainer,
480  AliAnalysisManager::GetCommonFileName());
481 
482  // Always connect input
483  mgr->ConnectInput(this, 0, mgr->GetCommonInputContainer());
484  mgr->ConnectOutput(this, 1, sumCon);
485 
486  return true;
487 }
488 //____________________________________________________________________
490 {
491  Printf("%s: %s", ClassName(), GetName());
492  Printf(" %22s: %d", "Scale by sin^2(theta)", fScaleDTheta);
493  Printf(" %22s: %f", "Delta phi shift", fDPhiShift);
494  Printf(" %22s: %f", "max Delta", fMaxDelta);
495  Printf(" %22s: %f", "Delta theta window", fDThetaWindow);
496  Printf(" %22s: %f", "Delta phi window", fDPhiWindow);
497  Printf(" %22s: %f", "phi overlap cut", fPhiOverlapCut);
498  Printf(" %22s: %f", "z-eta overlap cut", fZEtaOverlapCut);
499  Printf(" %22s: %s", "Filter strange",
500  fFilterMode==0 ? "no" :
501  fFilterMode==1 ? "random" : "track");
502 }
503 
504 //____________________________________________________________________
506 {
507  return *this;
508 }
509 //____________________________________________________________________
511 {
512  fContainer = new Container;
513  fContainer->SetOwner();
514  fContainer->SetName(GetName());
515  PostData(1, fContainer);
516 
517  // Create container of tracklets
518  if (WorkerInit()) return;
519 
520  AliWarning("Failed to initialize task on worker, making zombie");
521  SetZombie();
522 
523 }
524 
525 //____________________________________________________________________
527 {
528  AliInfo("Now initialising CDB");
529  AliAnalysisManager* anaMgr = AliAnalysisManager::GetAnalysisManager();
530  if (!anaMgr) {
531  AliError("No manager defined!");
532  return false;
533  }
534  // Check if we have the CDB connect task, and if so, do nothing as
535  // we rely on that task set up things properly.
536  const char* cdbNames[] = {"CDBconnect", "cdb", 0 };
537  const char** ptr = cdbNames;
538  while (*ptr) {
539  AliAnalysisTask* cdbConnect = anaMgr->GetTask(*ptr);
540  if (cdbConnect && cdbConnect->IsA()->InheritsFrom("AliTaskCDBconnect")) {
541  AliInfoF("CDB-connect task (%s: %s) present, do nothing",
542  cdbConnect->ClassName(), *ptr);
543  return true;
544  }
545  ptr++;
546  }
547  // Otherwise, we need to do stuff ourselves
548  AliInfo("Get the CDB manager");
549  AliCDBManager* cdbMgr = AliCDBManager::Instance();
550  if (!cdbMgr) {
551  AliError("Failed to get instance of CDB manager");
552  return false;
553  }
554  Int_t refRun = GetCDBReferenceRun();
555  TString refUrl = GetCDBReferenceURL();
556  AliWarningF("Using reference CDB storage \"%s\" and run \"%d\"",
557  refUrl.Data(), refRun);
558  // Set a very particular default storage. Perhaps we can do this
559  // with specific storages instead! Depends on whether the
560  // reconstruction also uses CDB - probably does - in which case we
561  // need specific storages for that too.
562  cdbMgr->SetDefaultStorage(refUrl);
563  // Now load our geometry - from a LHC10h run
564  AliInfo("Get Geometry entry");
565  AliCDBEntry* cdbEnt = cdbMgr->Get("GRP/Geometry/Data", refRun);
566  if (!cdbEnt) {
567  AliErrorF("No geometry found from %d", refRun);
568  return false;
569  }
570  // Initialize the geometry manager
571  AliInfo("Set Geometry");
572  AliGeomManager::SetGeometry(static_cast<TGeoManager*>(cdbEnt->GetObject()));
573  // Now perform mis-alignment - again based on an LHC10h run!
574  AliInfo("Misalign geometry");
575  if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",refRun,-1,-1)) {
576  AliErrorF("Failed to misalign geometry from %d", refRun);
577  return false;
578  }
579  return true;
580 }
581 //____________________________________________________________________
583 {
584  fTracklets = new TClonesArray(TrackletClassName());
585  // fTracklets->SetName(GetName());
586  TObject* obj = fTracklets;
587  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
588  AliAODHandler* ah =
589  dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
590  if (!ah) {
591  AliWarning("No AOD output handler set in analysis manager");
592  return false;
593  }
594  ah->AddBranch("TClonesArray", &obj);
595  // AliAODEvent* aod = ah->GetAOD();
596  // aod->Print();
597  // ah->GetTree()->Print();
598 
599  if (fFilterMode > 0) {
600  fStrangeLoss = new TParameter<double>("strLoss", 0);
601  ah->AddBranch("TParameter<double>", &fStrangeLoss);
602  }
603  return true;
604 }
605 
606 //____________________________________________________________________
608 {
609  fStatus = new TH1D("status","Processing status",kNStatus,0,kNStatus);
610  fStatus->SetDirectory(0);
611  fStatus->SetYTitle("Events");
612  fStatus->GetXaxis()->SetBinLabel(kSeen +1,"Seen");
613  fStatus->GetXaxis()->SetBinLabel(kClusters +1,"w/Clusters");
614  fStatus->GetXaxis()->SetBinLabel(kFilter +1,"Filtered");
615  fStatus->GetXaxis()->SetBinLabel(kIP +1,"w/IP");
616  fStatus->GetXaxis()->SetBinLabel(kReconstruct+1,"Reconstructed");
617  fStatus->GetXaxis()->SetBinLabel(kInject +1,"Injection");
618  fStatus->GetXaxis()->SetBinLabel(kTracklets +1,"Processed");
619  fStatus->GetXaxis()->SetBinLabel(kStored +1,"Stored");
620  fContainer->Add(fStatus);
621 
623  new TProfile2D("clustersVsTracklets",
624  "Correlation of clusters and tracklets",
625  1000,0,10000,
626  1000,0,10000);
627  fNClustersVsNTracklets->SetXTitle("#it{N}_{cluster,0}");
628  fNClustersVsNTracklets->SetYTitle("#it{N}_{cluster,1}");
629  fNClustersVsNTracklets->SetZTitle("#it{N}_{tracklets}");
630  fNClustersVsNTracklets->SetDirectory(0);
632 
633  if (!InitCDB()) return false;
634  if (!InitBranch()) return false;
635 
636  // PostData(1, fContainer);
637 
638  return true;
639 }
640 //____________________________________________________________________
642 {
643  // Clear container of tracklets
644  if (!fTracklets) {
645  AliWarning("Tracklets container not initialized, init must have failed!");
646  return;
647  }
648  fTracklets->Clear();
649 
650  if (!ProcessEvent()) return;
651 
652  MarkForStore();
653  PostData(1, fContainer);
654 }
655 
656 //____________________________________________________________________
658 {
659  AliVEvent* event = FindEvent();
660  if (!event) return false;
661  fStatus->Fill(kSeen);
662 
663  TTree* clusters = FindClusters();
664  const AliVVertex* ip = FindIP(event);
665  Bool_t ret = true;
666  if (fDebug > 0) AliInfoF("Cluster tree: %p, ip: %p", clusters, ip);
667  if (!ip) return false;
668  if (clusters) {
669  if (!HasGeometry()) ret = false;
670  if (fDebug > 1) AliInfoF("After geo check: %d", ret);
671  if (ret && !HasField(event)) ret = false;
672  if (fDebug > 1) AliInfoF("After field check: %d", ret);
673  // If we have clusters, then try to reconstruct the event (again)
674  if (ret) ret = Reconstruct(clusters, ip);
675  if (fDebug > 1) AliInfoF("After reconstruction: %d", ret);
676  }
677  else
678  // Otherwise, use the stored tracklets (max Delta = 1.5)
679  ret = ProcessTracklets(true, event->GetMultiplicity());
680 
681  CleanClusters(clusters);
682  if (fDebug > 1) AliInfoF("Return value: %d", ret);
683  return ret;
684 }
685 
686 //____________________________________________________________________
688 {
689  if (!AliGeomManager::GetGeometry()) {
690  AliError("No geometry loaded, needed for reconstruction");
691  AliError("Add the AliTaskCDBconnect to the train");
692  return false;
693  }
694  return true;
695 }
696 //____________________________________________________________________
698 {
699  if (!TGeoGlobalMagField::Instance()->GetField() &&
700  !event->InitMagneticField()) {
701  AliWarning("Failed to initialize magnetic field");
702  return false;
703  }
704  return true;
705 }
706 
707 //____________________________________________________________________
709 {
710  AliVEvent* event = InputEvent();
711  if (fDebug > 0 && !event) AliWarning("No event");
712  return event;
713 }
714 
715 //____________________________________________________________________
717 {
718  AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
719  AliVEventHandler* inh = mgr->GetInputEventHandler();
720  if (!inh->IsA()->InheritsFrom(AliESDInputHandlerRP::Class())) {
721  if (fDebug > 0)
722  AliWarningF("Clusters not available via input handler of class: %s",
723  inh->ClassName());
724  return 0;
725  }
726  AliESDInputHandlerRP* rph = static_cast<AliESDInputHandlerRP*>(inh);
727  TTree* tree = rph->GetTreeR("ITS");
728  if (!tree) {
729  AliWarning("Tree of clusters (rec.points) not found");
730  return 0;
731  }
732  fStatus->Fill(kClusters);
733  return FilterClusters(tree);
734 }
735 
736 //____________________________________________________________________
737 const AliVVertex* AliTrackletAODTask::FindIP(AliVEvent* event,
738  Double_t maxDispersion,
739  Double_t maxZError)
740 {
741  const AliVVertex* ip = event->GetPrimaryVertex();
742  if (!ip) return 0;
743  if (ip->GetNContributors() <= 0) {
744  if (fDebug > 0)
745  AliWarning("Not enough contributors for IP");
746  return 0;
747  }
748  // If this is from the Z vertexer, do some checks
749  if (ip->IsFromVertexerZ()) {
750  // Get covariance matrix
751  Double_t covar[6];
752  ip->GetCovarianceMatrix(covar);
753  Double_t sigmaZ = TMath::Sqrt(covar[5]);
754  if (sigmaZ >= maxZError) {
755  if (fDebug)
756  AliWarningF("IPz resolution = %f >= %f", sigmaZ, maxZError);
757  return 0;
758  }
759 
760  // If this IP doesn not derive from AliVertex, don't check dispersion.
761  if (ip->IsA()->InheritsFrom(AliVertex::Class())) {
762  const AliVertex* ipv = static_cast<const AliVertex*>(ip);
763  // Dispersion is the parameter used by the vertexer for finding the IP.
764  if (ipv->GetDispersion() >= maxDispersion) {
765  if (fDebug > 0)
766  AliWarningF("IP dispersion = %f >= %f",
767  ipv->GetDispersion(), maxDispersion);
768  return 0;
769  }
770  }
771  }
772 
773 #if 0
774  // If we get here, we either have a full 3D vertex or track
775  // vertex, and we should check if it is in range
776  if (ip->GetZ() < fIPzAxis.GetXmin() || ip->GetZ() > fIPzAxis.GetXmax()) {
777  AliWarningF("IPz = %fcm out of range [%f,%f]cm",
778  ip->GetZ(), fIPzAxis.GetXmin(), fIPzAxis.GetXmax());
779  return 0;
780  }
781 #endif
782  // Good vertex, return it
783  fStatus->Fill(kIP);
784  return ip;
785 }
786 
787 //____________________________________________________________________
789 {
790 
791  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
792  AliAODHandler* ah =
793  dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
794  if (!ah) return false;
795  ah->SetFillAOD(kTRUE);
796  if (fDebug > 0) AliInfo("Storing event");
797  fStatus->Fill(kStored);
798  return true;
799 }
800 
801 //____________________________________________________________________
802 Bool_t AliTrackletAODTask::Reconstruct(TTree* clusters, const AliVVertex* ip)
803 {
804  // Make an IP array for the reconstructor
805  Float_t ipv[3];
806  ipv[0] = ip->GetX();
807  ipv[1] = ip->GetY();
808  ipv[2] = ip->GetZ();
809  if (fDebug > 0)
810  AliInfoF("Reconstructing from cluster tree %p with (%4.2f,%4.2f,%7.4f",
811  clusters, ipv[0], ipv[1], ipv[2]);
812 
813  // --- Make the reconstructor --------------------------------------
814  fReco = 0;
815  AliITSMultRecBg reco;
816  reco.SetCreateClustersCopy (true);
817  reco.SetScaleDThetaBySin2T (fScaleDTheta);
818  reco.SetNStdDev (fMaxDelta);
819  reco.SetPhiWindow (fDPhiWindow);
820  reco.SetThetaWindow (fDThetaWindow);
821  reco.SetPhiShift (fDPhiShift);
822  reco.SetRemoveClustersFromOverlaps(fPhiOverlapCut > 0 ||
823  fZEtaOverlapCut > 0);
824  reco.SetPhiOverlapCut (fPhiOverlapCut);
825  reco.SetZetaOverlapCut (fZEtaOverlapCut);
826  reco.SetHistOn (false);
827  reco.SetRecType (AliITSMultRecBg::kData);
828 
829  // --- Run normal reconstruction -----------------------------------
830  reco.Run(clusters, ipv);
831  fStatus->Fill(kReconstruct);
832  fNClustersVsNTracklets->Fill(reco.GetNClustersLayer1(),
833  reco.GetNClustersLayer2(),
834  reco.GetMultiplicity()->GetNumberOfTracklets());
835 
836  // Save pointer to reco object - for MC combinatorics
837  fReco = &reco;
838  // And fill results into output branch
839  if (!ProcessTracklets(true, reco.GetMultiplicity())) {
840  AliWarning("Process tracklets (normal) failed");
841  return false;
842  }
843  fReco = 0;
844 
845  // --- Run again, but in injection mode ----------------------------
846  reco.SetRecType(AliITSMultRecBg::kBgInj);
847  reco.Run(clusters, ipv);
848  fStatus->Fill(kInject);
849 
850  // And fill results into output branch
851  if(!ProcessTracklets(false, reco.GetMultiplicity())) {
852  AliWarning("Process tracklets (injection) failed");
853  return false;
854  }
855  if (fDebug > 1) AliInfo("Returning true");
856  return true;
857 }
858 
859 //____________________________________________________________________
861  AliVMultiplicity* vmult)
862 {
863  // Now loop over all tracklets. Since this method is only called
864  // once per "reconstruction" pass, we save a lot of computing time
865  // since we loop over the tracklets of each reconstruction exactly
866  // once.
867  // if (!mult->IsA()->InheritsFrom(AliMultiplicity::Class())) return false;
868  AliMultiplicity* mult = static_cast<AliMultiplicity*>(vmult);
869  Int_t nTracklets = mult->GetNumberOfTracklets();
870  for (Int_t trackletNo = 0; trackletNo < nTracklets; trackletNo++) {
871  if (!ProcessTracklet(normal,mult,trackletNo)) return false;
872  }
873  if (normal) {
874  Int_t n0 = mult->GetNumberOfITSClusters(0);
875  Int_t n1 = mult->GetNumberOfITSClusters(1);
876  if (fDebug > 0)
877  Printf("%d x %d clusters -> %d", n0, n1, nTracklets);
878  if (n0 > 0 || n1 >> 0)
879  fNClustersVsNTracklets->Fill(n0, n1, nTracklets);
880  fStatus->Fill(kTracklets);
881  }
882  return true;
883 }
884 
885 //____________________________________________________________________
887 {
888  if (!fTracklets) return 0;
889  Int_t n = fTracklets->GetEntries();
890  return new((*fTracklets)[n]) AliAODTracklet;
891 }
892 //____________________________________________________________________
895  AliMultiplicity* mult,
896  Int_t no)
897 {
898  Double_t theta = mult->GetTheta (no);
899  Double_t phi = mult->GetPhi (no);
900  Double_t dTheta = mult->GetDeltaTheta(no);
901  Double_t dPhi = mult->GetDeltaPhi (no);
902  Double_t delta = mult->CalcDist (no);
903  AliAODTracklet* tracklet = MakeTracklet(normal);
904  tracklet->SetTheta (theta);
905  tracklet->SetPhi (phi);
906  tracklet->SetDTheta(dTheta);
907  tracklet->SetDPhi (dPhi);
908  tracklet->SetDelta (delta);
909  if (!normal) tracklet->SetInjection();
910 
911  return tracklet;
912 }
913 
914 /*********************************************************************
915  *
916  * Code for processing simulations
917  */
918 #ifndef __CINT__
919 #include <AliStack.h>
920 #include <AliMCEvent.h>
921 #include <AliGenEventHeader.h>
922 #include <TArrayF.h>
923 #include <TParticle.h>
924 #include <TParticlePDG.h>
925 #include <TDatabasePDG.h>
926 #include <AliITSgeomTGeo.h>
927 #include <TRandom.h>
928 #include <TROOT.h>
929 #else
930 class TParticlePDG;
931 class TParticle;
932 class AliStack;
933 #endif
934 //====================================================================
941 {
942 public:
947  : AliTrackletAODTask(),
948  fFilterWeights(0),
949  fSeenTrackPDGs(0),
950  fUsedTrackPDGs(0),
951  fSeenClusterPDGs(0),
953  {}
959  AliTrackletAODMCTask(const char* name)
960  : AliTrackletAODTask(name),
961  fFilterWeights(0),
962  fSeenTrackPDGs(0),
963  fUsedTrackPDGs(0),
964  fSeenClusterPDGs(0),
966  {}
973  : AliTrackletAODTask(other),
974  fFilterWeights(0),
975  fSeenTrackPDGs(0),
976  fUsedTrackPDGs(0),
977  fSeenClusterPDGs(0),
979  {}
988  {
989  return *this;
990  }
997 protected:
1002  /*
1003  * Initialize the worker
1004  */
1005  Bool_t WorkerInit();
1012  virtual Bool_t ProcessEvent();
1019  /* @} */
1029  virtual TTree* FilterClusters(TTree* t);
1042  virtual void FilterClustersRandom(TTree* t, TTree* copy,
1043  Double_t cent, Double_t ipz);
1057  virtual void FilterClustersTrack(TTree* t, TTree* copy,
1058  Double_t cent, Double_t ipz);
1064  virtual void CleanClusters(TTree*& t);
1072  Double_t LookupWeight(TParticle* particle,
1073  Double_t cent=0,
1074  Double_t ipz=0) const;
1082  Bool_t KeepIt(Double_t weight) const;
1083  /* @} */
1093  virtual const char* GetCDBReferenceURL() const
1094  {
1095  return "alien://Folder=/alice/simulation/2008/v4-15-Release/Residual";
1096  }
1097  /* @} */
1107  virtual const char* TrackletClassName() const { return "AliAODMCTracklet"; }
1115  virtual AliAODTracklet* MakeTracklet(Bool_t normal);
1125  virtual AliAODTracklet* ProcessTracklet(Bool_t normal,
1126  AliMultiplicity* mult,
1127  Int_t no);
1128  /* @} */
1142  Int_t FindPrimaryParentID(Int_t label) const;
1151  TParticle* FindPrimaryParent(Int_t label) const;
1160  Int_t FindParent(Int_t label) const;
1175  Int_t FindParents(Int_t label, TArrayI& fill, Int_t offset) const;
1186  Int_t CommonParent(Int_t label, const TArrayI& fill) const;
1187  /* @} */
1192 
1194 
1196 
1198 
1199  ClassDef(AliTrackletAODMCTask,1);
1200 };
1201 
1202 #define ADDP(P,M,L,Q) \
1203  db->AddParticle("p" #P, "p" #P,M,false,L,Q,"Baryon",P); \
1204  db->AddAntiParticle("p" #P "_bar", -P)
1205 
1206 //____________________________________________________________________
1208 {
1209  if (!AliTrackletAODTask::WorkerInit()) return false;
1210 
1211  const TAxis& pdgAxis = PdgAxis();
1212  TAxis layerAxis(2,0,2);
1213  layerAxis.SetBinLabel(1,"Layer 0");
1214  layerAxis.SetBinLabel(2,"Layer 1");
1215  FixAxis(layerAxis,"Layer");
1216  fSeenTrackPDGs = Make1D(fContainer, "seenTrackPdg",
1217  "Seen track PDGs", kGreen+1,20, pdgAxis);
1218  fUsedTrackPDGs = Make1D(fContainer, "usedTrackPdg",
1219  "Used track PDGs", kBlue+1,24, pdgAxis);
1220  fSeenClusterPDGs = Make2D(fContainer, "seenClusterPdg",
1221  "Seen cluster PDGs", kGreen+1,21,pdgAxis,layerAxis);
1222  fUsedClusterPDGs = Make2D(fContainer, "usedClusterPdg",
1223  "Used cluster PDGs", kBlue+1,25,pdgAxis,layerAxis);
1224 
1225 
1226  // Register additional particles from EPOS-LHC
1227  TDatabasePDG* db = TDatabasePDG::Instance();
1228  db->GetParticle(2212);
1229  db->AddParticle("deut","deut",1.876560,false,0.000000,3,"Baryon",1000010020);
1230  db->AddAntiParticle("deut_bar",-1000010020);
1231  db->AddParticle("trit","trit",2.816700,false,0.000000,3,"Baryon",1000010030);
1232  db->AddAntiParticle("trit_bar",-1000010030);
1233  db->AddParticle("alph","alph",3.755000,false,0.000000,6,"Baryon",1000020040);
1234  db->AddAntiParticle("alph_bar",-1000020040);
1235  ADDP(32224, 1.524000,0.145000,6);
1236  ADDP(12224, 1.816000,0.350000,6);
1237  ADDP(12222, 2.108000,0.350000,6);
1238  ADDP(12212, 1.525720,0.145000,3);
1239  ADDP(2124, 1.819440,0.150000,3);
1240  ADDP(32214, 2.113160,0.150000,3);
1241  ADDP(32212, 2.406880,0.145000,3);
1242  ADDP(12214, 2.700600,0.300000,3);
1243  ADDP(22124, 2.994320,0.150000,3);
1244  ADDP(12122, 3.288040,0.300000,3);
1245  ADDP(13222, 1.575200,0.080000,3);
1246  ADDP(23222, 1.768100,0.100000,3);
1247  ADDP(13226, 1.961000,0.170000,3);
1248  ADDP(12112, 1.524430,0.350000,0);
1249  ADDP(1214, 1.816860,0.150000,0);
1250  ADDP(32114, 2.109290,0.150000,0);
1251  ADDP(32112, 2.401720,0.145000,0);
1252  ADDP(12114, 2.694150,0.300000,0);
1253  ADDP(21214, 2.986579,0.150000,0);
1254  ADDP(11212, 3.279010,0.300000,0);
1255  ADDP(13122, 1.761000,0.040000,0);
1256  ADDP(3124, 1.950500,0.016000,0);
1257  ADDP(23122, 2.140000,0.090000,0);
1258  ADDP(13212, 2.329500,0.080000,0);
1259  ADDP(23212, 2.519000,0.100000,0);
1260  ADDP(43122, 2.708500,0.145000,0);
1261  ADDP(13216, 2.898000,0.170000,0);
1262  ADDP(31114, 1.524000,0.145000,-3);
1263  ADDP(11114, 1.816000,0.350000,-3);
1264  ADDP(11112, 2.108000,0.350000,-3);
1265  ADDP(13112, 1.577600,0.080000,-3);
1266  ADDP(23112, 1.767700,0.100000,-3);
1267  ADDP(13116, 1.957800,0.170000,-3);
1268  ADDP(9900110,0.000000,0.000000,0);
1269  ADDP(9900210,0.000000,0.000000,0);
1270  ADDP(9900220,0.000000,0.000000,0);
1271  ADDP(9900330,0.000000,0.000000,0);
1272  ADDP(9900440,0.000000,0.000000,0);
1273  ADDP(9902210,0.000000,0.000000,0);
1274  ADDP(9902110,0.000000,0.000000,0);
1275  ADDP(88, 0.000000,0.000000,0);
1276  ADDP(90, 0.000000,0.000000,0);
1277  ADDP(990, 0.000000,0.000000,0);
1278  ADDP(99999, 0.000000,0.000000,0);
1279 
1280  return true;
1281 }
1282 
1283 //____________________________________________________________________
1285 {
1286  if (fStrangeLoss) fStrangeLoss->SetVal(0);
1287  if (!t || fFilterMode <= 0) {
1288  if (fDebug > 0) AliInfo("Returning original cluster tree");
1289  return t;
1290  }
1291 
1292  TDirectory* savDir = gDirectory;
1293  gDirectory = gROOT;
1294  TTree* copy = new TTree("TreeR", "TreeR");
1295  copy->SetAutoFlush(0);// Keep in memory
1296 
1297  Double_t cent = 0;
1298  Double_t ipz = 0;
1299  if (MCEvent()->GenEventHeader()) {
1300  TArrayF v(3);
1301  MCEvent()->GenEventHeader()->PrimaryVertex(v);
1302  ipz = v[2];
1303  }
1304 
1305  if (fFilterMode == 1) FilterClustersRandom(t, copy, cent, ipz);
1306  else FilterClustersTrack(t, copy, cent, ipz);
1307 
1308  savDir->cd();
1309  if (fDebug > 0)
1310  AliInfoF("Returning reduced cluster tree %p (original %p)", copy, t);
1311  fStatus->Fill(kFilter);
1312  return copy;
1313 }
1314 
1315 
1316 //____________________________________________________________________
1318  TTree* copy,
1319  Double_t cent,
1320  Double_t ipz)
1321 {
1322  TClonesArray* in = new TClonesArray("AliITSRecPoint");
1323  TClonesArray* out = new TClonesArray("AliITSRecPoint");
1324  Int_t outN = 0;
1325  copy->Branch("ITSRecPoints", &out);
1326  t->SetBranchAddress("ITSRecPoints", &in);
1327 
1328  // Printf("Filtering clusters from strange particles");
1329 
1330  Int_t min1 = AliITSgeomTGeo::GetModuleIndex(1,1,1);
1331  Int_t max1 = AliITSgeomTGeo::GetModuleIndex(2,1,1);
1332  Int_t min2 = AliITSgeomTGeo::GetModuleIndex(2,1,1);
1333  Int_t max2 = AliITSgeomTGeo::GetModuleIndex(3,1,1);
1334  Int_t nTotal = 0;
1335  Int_t nKept = 0;
1336 
1337  // Loop over the modules of the SPD
1338  for (Int_t i = 0; i < max2; i++) {
1339  in->Clear();
1340  out->Clear();
1341  outN = 0;
1342 
1343  // Read in module data
1344  t->GetEntry(i);
1345 
1346  // Loop over all clusters in the module
1347  Int_t inN = in->GetEntries();
1348  for (Int_t j = 0; j < inN; j++) {
1349  AliITSRecPoint* inCl = static_cast<AliITSRecPoint*>(in->At(j));
1350  if (!inCl) continue;
1351 
1352  // Loop over labels in the module
1353  Double_t weight = 1;
1354  for (Int_t k = 0; k < 3; k++) {
1355  Int_t label = inCl->GetLabel(k);
1356  if (label <= 0) continue;
1357 
1358  // Check primary parent particle type
1359  TParticle* parent = FindPrimaryParent(label);
1360  if (!parent) continue;
1361 
1362  // weight *= GetWeight(parent->GetPdgCode());
1363  weight *= LookupWeight(parent,cent,ipz);
1364  if (weight > 1)
1365  Printf("Cluster from %+5d: %7.5f", parent->GetPdgCode(), weight);
1366  // Printf("Parent %d is a K^0_S: %d", k+1, parent->GetPdgCode());
1367  // Randomly remove clusters from the right kind of primary
1368  // parent particles
1369  }
1370  if (weight > 1) {
1371  nTotal++;
1372  if (!KeepIt(weight)) continue;
1373  nKept++;
1374  }
1375  new ((*out)[outN++]) AliITSRecPoint(*inCl);
1376  }
1377  // Printf("Kept %d out of %d clusters", outN, inN);
1378  copy->Fill();
1379  }
1380  if (nTotal > 0) {
1381  Double_t loss = 1-Float_t(nKept)/nTotal;
1382  Printf("Kept %d out of %d clusters from strange primaries (%4.1f%%)",
1383  nKept, nTotal, 100.*loss);
1384  if (fStrangeLoss) fStrangeLoss->SetVal(loss);
1385  }
1386  delete in;
1387  delete out;
1388 }
1389 
1390 //____________________________________________________________________
1392  TTree* copy,
1393  Double_t cent,
1394  Double_t ipz)
1395 {
1396  TClonesArray* in = new TClonesArray("AliITSRecPoint");
1397  TClonesArray* out = new TClonesArray("AliITSRecPoint");
1398  Int_t outN = 0;
1399  copy->Branch("ITSRecPoints", &out);
1400  t->SetBranchAddress("ITSRecPoints", &in);
1401 
1402  // Loop over the stack
1403  Int_t nTotal = 0;
1404  Int_t nKept = 0;
1405  AliStack* stack = MCEvent()->Stack();
1406  Int_t nTracks = stack->GetNtrack();
1407  TBits kept(nTracks);
1408  TBits seen(nTracks);
1409  kept.ResetAllBits(true);
1410  seen.ResetAllBits(false);
1411 
1412  for (Int_t trackNo = nTracks; trackNo--; ) {
1413  // Get the primary parent identifier
1414  Int_t parent = FindPrimaryParentID(trackNo);
1415  if (parent < 0) continue;
1416  // Check if we've seen this parent already
1417  if (seen.TestBitNumber(parent)) continue;
1418  seen.SetBitNumber(parent,true);
1419 
1420  // Get the primary parent track and weight
1421  TParticle* par = stack->Particle(parent);
1422  // Double_t weight = GetWeight(par->GetPdgCode());
1423  Double_t weight = LookupWeight(par,cent,ipz);
1424  Bool_t keep = true;
1425  if (weight > 1) {
1426  nTotal++;
1427  keep = KeepIt(weight);
1428  if (keep) nKept++;
1429  }
1430  Int_t bin = PdgBin(par->GetPdgCode());
1431  fSeenTrackPDGs->Fill(bin);
1432  if (keep) fUsedTrackPDGs->Fill(bin);
1433  if (fDebug > 1)
1434  Printf("Primary parent %6d from a %6d %6s (%f)",
1435  parent, par->GetPdgCode(), keep ? "kept" : "marked", weight);
1436  kept.SetBitNumber(parent, keep);
1437  par->SetWeight(keep ? 1 : weight);
1438  }
1439  // At this point, we have investigated all relevant primary
1440  // particles and thrown a dice to see if we should remove clusters
1441  // that correspond to these primary particles.
1442  if (fDebug > 0 && nTotal > 0) {
1443  Double_t loss = 1-Float_t(nKept)/nTotal;
1444  Printf("Kept %d out of %d strange primaries (%4.1f%% loss)",
1445  nKept, nTotal, 100.*loss);
1446  if (fStrangeLoss) fStrangeLoss->SetVal(loss);
1447  }
1448 
1449  // Next thing is to actually remove the clusters
1450  // Printf("Filtering clusters from strange particles");
1451  Int_t min1 = AliITSgeomTGeo::GetModuleIndex(1,1,1);
1452  Int_t max1 = AliITSgeomTGeo::GetModuleIndex(2,1,1);
1453  Int_t min2 = AliITSgeomTGeo::GetModuleIndex(2,1,1);
1454  Int_t max2 = AliITSgeomTGeo::GetModuleIndex(3,1,1);
1455  Int_t inT = 0;
1456  Int_t outT = 0;
1457  // Loop over the modules of the SPD
1458  for (Int_t i = 0; i < max2; i++) {
1459  in->Clear();
1460  out->Clear();
1461  outN = 0;
1462 
1463  // Read in module data
1464  t->GetEntry(i);
1465 
1466  // Loop over all clusters in the module
1467  Int_t inN = in->GetEntries();
1468  inT += inN;
1469  for (Int_t j = 0; j < inN; j++) {
1470  AliITSRecPoint* inCl = static_cast<AliITSRecPoint*>(in->At(j));
1471  if (!inCl) continue;
1472 
1473  // Loop over labels of the cluster
1474  Bool_t toRemove = false;
1475  Int_t pdg = 0;
1476  for (Int_t k = 0; k < 3; k++) {
1477  Int_t label = inCl->GetLabel(k);
1478  if (label <= 0) continue;
1479 
1480  // Check primary parent particle type
1481  Int_t parent = FindPrimaryParentID(label);
1482  if (parent < 0) continue;
1483 
1484  if (pdg == 0) pdg = stack->Particle(parent)->GetPdgCode();
1485  if (!kept.TestBitNumber(parent)) toRemove = true;
1486  if (fDebug > 3) {
1487  Printf("Cluster %6d from parent %6d (%7.5f) of type %6d %s",
1488  j, parent,
1489  stack->Particle(parent)->GetPdgCode(),
1490  stack->Particle(parent)->GetWeight(),
1491  toRemove ? "removed" : "kept");
1492  }
1493  }
1494  // We've now looked at all primaries for this cluster, and if
1495  // any of these primaries was flagged for removal, then we
1496  // should remove the cluster.
1497  if (fDebug > 3)
1498  Printf("Cluster %2d/%3d from %4d is %s",
1499  i, j, pdg, toRemove ? "removed" : "kept");
1500  fSeenClusterPDGs->Fill(PdgBin(pdg), i < max1);
1501  if (toRemove) continue;
1502  fUsedClusterPDGs->Fill(PdgBin(pdg), i < max1);
1503  new ((*out)[outN++]) AliITSRecPoint(*inCl);
1504  }
1505  // Printf("Kept %d out of %d clusters", outN, inN);
1506  outT += outN;
1507  copy->Fill();
1508  }
1509  if (fDebug > 0)
1510  Printf("Wrote out %6d out of %6d clusters (%4.1f%% loss)",
1511  outT, inT, (inT > 0 ? 100*float(inT-outT)/inT : 100));
1512  delete in;
1513  delete out;
1514 }
1515 
1516 //____________________________________________________________________
1518  Double_t cent,
1519  Double_t ipz) const
1520 {
1521 #if 1
1522  if (!fFilterWeights) return 1;
1523  return fFilterWeights->LookupWeight(particle,cent,ipz);
1524 #else
1525  const Double_t k0s = 1.52233299626516083e+00; // 310 - K^0_S weight
1526  const Double_t kpm = (1.43744204476109627e+00*
1527  9.82150320171071400e-01); // 321 - K^{+/-}
1528  const Double_t lam = 2.75002089647900005e+00; // 3122 - lambda
1529  const Double_t sig = 2.75002089647899961e+00; // 3212 - sigma
1530  const Double_t xi = 3.24109605656453548e+00; // 3322 - Xi
1531 
1532  switch (TMath::Abs(particle->GetPdgCode())) {
1533  case 310: return k0s;
1534  case 321: return kpm;
1535  case 3122: return lam;
1536  // case 3212: return sig; // Old, wrong code
1537  case 3112:
1538  case 3222: return sig;
1539  case 3312: return xi;
1540  // case 3322: return xi; // Old, wrong code
1541  }
1542  return 1;
1543  }
1544 #endif
1545 }
1546 
1547 //____________________________________________________________________
1549 {
1550  if (weight <= 1) return true;
1551  Double_t chance = 1 - 1 / weight; // chance is 1 minus inverse
1552  return (gRandom->Uniform() >= chance);
1553 }
1554 
1555 //____________________________________________________________________
1557 {
1558  if (!t || fFilterMode <= 0) return;
1559 
1560  delete t;
1561  t = 0;
1562 }
1563 
1564 //____________________________________________________________________
1566 {
1567  if (!AliTrackletAODTask::ProcessEvent()) return false;
1568 
1569  // create generated "tracklets"
1570  if (!ProcessGenerated()) return false;
1571 
1572  return true;
1573 }
1574 
1575 //____________________________________________________________________
1577 {
1578  if (fDebug > 0) AliInfo("Processing generated particles");
1579  AliStack* stack = MCEvent()->Stack();
1580  Int_t nTracks = stack->GetNtrack();
1581  for (Int_t trackNo = nTracks; trackNo--; ) {
1582  if (!stack->IsPhysicalPrimary(trackNo)) {
1583  // AliWarningF("Track # %6d not a primary", trackNo);
1584  // Not a primary, go on
1585  continue;
1586  }
1587  // Get the particle
1588  TParticle* particle = stack->Particle(trackNo);
1589  if (!particle) {
1590  AliWarningF("No particle found for track # %d", trackNo);
1591  continue;
1592  }
1593  TParticlePDG* pdg = particle->GetPDG();
1594  if (!pdg) {
1595  AliWarningF("Unknown PDG code: %d", particle->GetPdgCode());
1596  // continue; // Unknown particle
1597  }
1598  // if (pdg->Charge() == 0) {
1599  // Uncharged - don't care
1600  // continue;
1601  // }
1602 
1603  // Get theta
1604  Double_t theta = particle->Theta();
1605  // Check for beam-like particle
1606  if (theta < 1e-6 || TMath::Abs(theta-TMath::Pi()) < 1e-6) {
1607  if (fDebug > 0)
1608  AliWarningF("Track # %6d is beam-like (%f)", trackNo,
1609  TMath::RadToDeg()*theta);
1610  continue;
1611  }
1612  Double_t eta = -TMath::Log(TMath::ATan(theta/2));
1613  // If the pseudorapidity is way beyond the SPD acceptance, do not
1614  // write down this generated "tracklet" - to save space.
1615  if (TMath::Abs(eta) > 3) continue;
1616 
1617  Double_t phi = particle->Phi();
1618  AliAODTracklet* tracklet = MakeTracklet(false);
1619  AliAODMCTracklet* mc = static_cast<AliAODMCTracklet*>(tracklet);
1620  mc->SetGenerated();
1621  mc->SetTheta(theta);
1622  mc->SetPhi(phi);
1623  mc->SetDTheta(0);
1624  mc->SetDPhi (0);
1625  mc->SetDelta (0);
1626  mc->SetParentPdg(particle->GetPdgCode());
1627  mc->SetParentPt (particle->Pt());
1628  if (particle->GetWeight() > 1) mc->SetSuppressed();
1629  if (pdg && pdg->Charge() == 0) mc->SetNeutral();
1630  }
1631  if (fDebug > 2) AliInfo("Returning true from generated");
1632  return true;
1633 }
1634 
1635 //____________________________________________________________________
1637 {
1638  if (!fTracklets) return 0;
1639  Int_t n = fTracklets->GetEntries();
1640  return new((*fTracklets)[n]) AliAODMCTracklet;
1641 }
1642 //____________________________________________________________________
1645  AliMultiplicity* mult,
1646  Int_t no)
1647 {
1648  AliAODTracklet* tracklet =
1649  AliTrackletAODTask::ProcessTracklet(normal,mult,no);
1650  if (!normal) return tracklet;
1651 
1652  AliAODMCTracklet* mc = static_cast<AliAODMCTracklet*>(tracklet);
1653  Int_t label0 = mult->GetLabel(no, 0);
1654  Int_t label1 = mult->GetLabel(no, 1);
1655  TParticle* parent0 = FindPrimaryParent(label0);
1656  if (parent0) {
1657  mc->SetParentPdg(parent0->GetPdgCode());
1658  mc->SetParentPt (parent0->Pt());
1659  }
1660  if (label0 != label1) {
1661  TParticle* parent1 = FindPrimaryParent(label1);
1662  if (parent1) {
1663  mc->SetParentPdg(parent1->GetPdgCode(), true);
1664  mc->SetParentPt (parent1->Pt(), true);
1665  }
1666  mc->SetCombinatorics();
1667  // Here, we could track back in the cluster labels to see if we
1668  // have the same ultimate mother of both clusters. Since this
1669  // isn't really used, we do not do that
1670  if (fReco) {
1671  TArrayI parents(50); // At most 50 levels deep
1672  // Tracklet parameters
1673  Float_t* ftrack = fReco->GetTracklet(no);
1674  // Cluster identifiers
1675  Int_t clus0 = Int_t(ftrack[AliITSMultReconstructor::kClID1]);
1676  Int_t clus1 = Int_t(ftrack[AliITSMultReconstructor::kClID2]);
1677  // Cluster labelsx
1678  Float_t* fclus0 = (fReco->GetClusterOfLayer(0,clus0) +
1679  AliITSMultReconstructor::kClMC0);
1680  Float_t* fclus1 = (fReco->GetClusterOfLayer(1,clus1) +
1681  AliITSMultReconstructor::kClMC0);
1682  // Loop over three inner layer cluster labels
1683  Int_t offset = 0;
1684  for (Int_t i = 0; i < 3; i++)
1685  offset = FindParents(Int_t(fclus0[i]), parents, offset);
1686  // Loop over three outer layer cluster labels
1687  Bool_t distinct = true;
1688  for (Int_t i = 0; i < 3; i++) {
1689  Float_t flbl = fclus1[i];
1690  // Be careful not to get overflows on conversions
1691  if (flbl > nextafter(INT_MAX, 0) || flbl < nextafter(INT_MIN, 0))
1692  continue; // Would overflow
1693  if (CommonParent(Int_t(flbl), parents)) {
1694  distinct = false;
1695  // We break out as soon as we find a common parent.
1696  break;
1697  }
1698  } // loop over outer layer
1699  if (distinct) mc->SetDistinct();
1700  } // if (fReco)
1701  }
1702  else {
1703  if (!MCEvent()->Stack()->IsPhysicalPrimary(label0))
1704  mc->SetSecondary();
1705  }
1706  return tracklet;
1707 }
1708 
1709 //____________________________________________________________________
1711 {
1712  AliStack* stack = MCEvent()->Stack();
1713  Int_t nTracks = stack->GetNtrack();
1714  TParticle* particle = stack->Particle(label);
1715  if (!particle) return -1;
1716  Int_t ret = particle->GetFirstMother();
1717  if (ret > nTracks) return -1;
1718  return ret;
1719 }
1720 
1721 //____________________________________________________________________
1723  TArrayI& fill,
1724  Int_t offset) const
1725 {
1726  // If offset is negative, then we don't really store the labels, but
1727  // only check against those already stored.
1728  Int_t i = offset;
1729  Int_t lbl = label;
1730  while (i < fill.GetSize()-1 && lbl >= 1) {
1731  fill[i] = lbl;
1732  i++;
1733  lbl = FindParent(lbl);
1734  }
1735  // If we get here, and we're checking, that means we did not find a
1736  // common ancestor, so we return 0 (false). If we're not checking,
1737  // then we return the number of elements assigned in the passed
1738  // cache.
1739  return i;
1740 }
1741 
1742 //____________________________________________________________________
1744  const TArrayI& fill) const
1745 {
1746  Int_t i = 0;
1747  Int_t lbl = label;
1748  while (i < fill.GetSize()-1 && lbl >= 1) {
1749  // If we're checking, just see if we have the label in the list
1750  // already. If we do, then return the index
1751  for (Int_t j = 0; j < fill.GetSize(); j++) {
1752  if (fill[j] == lbl) return j;
1753  }
1754  i++;
1755  lbl = FindParent(lbl);
1756  }
1757  // If we get here, and we're checking, that means we did not find a
1758  // common ancestor, so we return 0 (false). If we're not checking,
1759  // then we return the number of elements assigned in the passed
1760  // cache.
1761  return false;
1762 }
1763 
1764 
1765 //____________________________________________________________________
1767 {
1768  Int_t parent = FindPrimaryParentID(label);
1769  if (parent < 0) return 0;
1770  return MCEvent()->Stack()->Particle(parent);
1771 }
1772 //____________________________________________________________________
1774 {
1775  AliStack* stack = MCEvent()->Stack();
1776  Int_t nTracks = stack->GetNtrack();
1777  Int_t trackNo = label;
1778  if (trackNo > nTracks || trackNo <= 0) return 0;
1779  TParticle* particle = stack->Particle(label);
1780  while (!stack->IsPhysicalPrimary(trackNo)) {
1781  trackNo = particle->GetFirstMother();
1782  // If we have hit the top
1783  if (trackNo < 0) return 0;
1784  // Partice first next iteration
1785  particle = stack->Particle(trackNo);
1786  }
1787  return trackNo;
1788 }
1789 //====================================================================
1791 {
1792  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1793  if (!mgr) {
1794  ::Error("Create","No analysis manager to connect to.");
1795  return 0;
1796  }
1797  AliAODHandler* ah =
1798  dynamic_cast<AliAODHandler*>(mgr->GetOutputEventHandler());
1799  if (!ah) {
1800  ::Error("Create","No AOD output handler!");
1801  return 0;
1802  }
1803 
1804  Bool_t mc = mgr->GetMCtruthEventHandler() != 0;
1805  AliTrackletAODTask* ret = 0;
1806  if (mc) ret = new AliTrackletAODMCTask("MidRapidityMC");
1807  else ret = new AliTrackletAODTask("MidRapidity");
1808  if (ret) ret->Connect();
1809 
1810  if (weights && weights[0] != '\0') {
1811  TUrl wurl(weights);
1812  TFile* wfile = TFile::Open(wurl.GetFile());
1813  if (!wfile) {
1814  ::Warning("Create", "Failed to open weights file: %s",
1815  wurl.GetUrl());
1816  return 0;
1817  }
1818  TString wnam(wurl.GetAnchor());
1819  if (wnam.IsNull()) wnam = "weights";
1820 
1821  TObject* wobj = wfile->Get(wnam);
1822  if (!wobj) {
1823  ::Warning("Create", "Failed to get weights %s from file %s",
1824  wnam.Data(), wfile->GetName());
1825  return 0;
1826  }
1827  if (!wobj->IsA()->InheritsFrom(AliTrackletBaseWeights::Class())) {
1828  ::Warning("Create", "Object %s from file %s not an "
1829  "AliTrackletBaseWeights but a %s",
1830  wnam.Data(), wfile->GetName(), wobj->ClassName());
1831  return 0;
1832  }
1833  ret->SetFilterWeights(static_cast<AliTrackletBaseWeights*>(wobj));
1834  }
1835  return ret;
1836 }
1837 
1838 
1839 //====================================================================
1840 
1841 //
1842 // EOF
1843 //
1844 
1845 
Int_t pdg
static const TAxis & PdgAxis()
void SetPhiOverlapCut(Double_t x=0.005)
void SetPhi(Real_t x)
static TH2 * Make2D(Container *c, const TString &name, const TString &title, Color_t color, Style_t style, const TAxis &xAxis, const TAxis &yAxis)
double Double_t
Definition: External.C:58
void UserExec(Option_t *)
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
static AliTrackletAODTask * Create(const char *weights="")
void SetParentPdg(Short_t pdg, Bool_t second=false)
AliTrackletAODTask & operator=(const AliTrackletAODTask &other)
Bool_t KeepIt(Double_t weight) const
void SetTheta(Real_t x)
AliTrackletAODMCTask(const char *name)
AliITSMultRecBg * fReco
void SetDPhiWindow(Double_t x=0.06)
AliStack * stack
void SetDTheta(Real_t x)
TRandom * gRandom
void SetDPhiShift(Double_t x=0.0045)
virtual void CleanClusters(TTree *&t)
Definition: External.C:92
virtual Bool_t ProcessEvent()
TString kData
Declare data MC or deltaAOD.
virtual void SetFilterWeights(AliTrackletBaseWeights *w)
Utilities for midrapidity analysis.
Tracklet AOD object.
virtual TTree * FilterClusters(TTree *t)
TParticle * FindPrimaryParent(Int_t label) const
virtual void SetFilterWeights(AliTrackletBaseWeights *w)
virtual const char * TrackletClassName() const
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)
Encode simulation weights for 2nd pass.
Definition: External.C:212
Bool_t ProcessTracklets(Bool_t normal, AliVMultiplicity *mult)
virtual void FilterClustersTrack(TTree *t, TTree *copy, Double_t cent, Double_t ipz)
TParameter< double > * fStrangeLoss
virtual Bool_t ProcessEvent()
Int_t mode
Definition: anaM.C:41
virtual Int_t GetCDBReferenceRun() const
void SetParentPt(Real_t pt, Bool_t second=false)
Double_t LookupWeight(AliAODTracklet *tracklet, Double_t cent, Double_t ipz, TH2 *corr=0) const
AliTrackletBaseWeights * fFilterWeights
virtual AliAODTracklet * MakeTracklet(Bool_t normal)
TProfile2D * fNClustersVsNTracklets
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
Definition: External.C:220
void SetMaxDelta(Double_t x=25)
virtual AliAODTracklet * ProcessTracklet(Bool_t normal, AliMultiplicity *mult, Int_t no)
virtual TTree * FilterClusters(TTree *t)
static Int_t PdgBin(Int_t pdg)
void Print(Option_t *) const
virtual const char * GetCDBReferenceURL() const
AliTrackletAODMCTask(const AliTrackletAODMCTask &other)
const char Option_t
Definition: External.C:48
virtual void SetFilterMode(Int_t mode)
Bool_t Reconstruct(TTree *clusters, const AliVVertex *ip)
AliTrackletAODMCTask & operator=(const AliTrackletAODMCTask &other)
bool Bool_t
Definition: External.C:53
Double_t LookupWeight(TParticle *particle, Double_t cent=0, Double_t ipz=0) const
static TH1 * Make1D(Container *c, const TString &name, const TString &title, Color_t color, Style_t style, const TAxis &xAxis)
void SetDelta(Real_t x)
void SetDPhi(Real_t x)
static void FixAxis(TAxis &axis, const char *title=0)
Int_t FindParents(Int_t label, TArrayI &fill, Int_t offset) const
Definition: External.C:196
#define ADDP(P, M, L, Q)
virtual void FilterClustersRandom(TTree *t, TTree *copy, Double_t cent, Double_t ipz)
virtual const char * TrackletClassName() const
void Terminate(Option_t *)