AliPhysics  e34b7ac (e34b7ac)
 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; }
179  virtual void SetFilterK0S(Bool_t filter) { fFilterK0S = filter; }
180  /* @} */
181 protected:
189  virtual Bool_t WorkerInit();
195  Bool_t InitCDB();
201  virtual Int_t GetCDBReferenceRun() const { return 137161; }
207  virtual const char* GetCDBReferenceURL() const
208  {
209  return "alien://Folder=/alice/data/2010/OCDB";
210  }
216  Bool_t InitBranch();
227  virtual Bool_t ProcessEvent();
239  Bool_t HasField(AliVEvent* event);
245  TTree* FindClusters();
251  virtual TTree* FilterClusters(TTree* t) { return t; }
257  virtual void CleanClusters(TTree*& t) {};
263  AliVEvent* FindEvent();
273  const AliVVertex* FindIP(AliVEvent* event,
274  Double_t maxDispersion=0.04,
275  Double_t maxZError=0.25);
282  Bool_t Reconstruct(TTree* clusters, const AliVVertex* ip);
289  /* @} */
299  virtual const char* TrackletClassName() const { return "AliAODTracklet"; }
307  virtual AliAODTracklet* MakeTracklet(Bool_t normal);
316  Bool_t ProcessTracklets(Bool_t normal, AliVMultiplicity* mult);
326  virtual AliAODTracklet* ProcessTracklet(Bool_t normal,
327  AliMultiplicity* mult,
328  Int_t no);
329  /* @} */
330 
332  TClonesArray* fTracklets;
333  /* Output container */
334  // Container* fContainer;
350  AliITSMultRecBg* fReco;
351 
355 
357 };
358 //____________________________________________________________________
360  : AliAnalysisTaskSE(),
361  fTracklets(0),
362  // fContainer(0),
363  fScaleDTheta(true),
364  fMaxDelta(25),
365  fDPhiShift(0.0045),
366  fDThetaWindow(0.025),
367  fDPhiWindow(0.06),
368  fPhiOverlapCut(0.005),
369  fZEtaOverlapCut(0.05),
370  fReco(0),
371  fFilterK0S(false),
372  fK0SLoss(0)
373 {}
374 //____________________________________________________________________
376  : AliAnalysisTaskSE("tracklets"),
377  fTracklets(0),
378  // fContainer(0),
379  fScaleDTheta(true),
380  fMaxDelta(25),
381  fDPhiShift(0.0045),
382  fDThetaWindow(0.025),
383  fDPhiWindow(0.06),
384  fPhiOverlapCut(0.005),
385  fZEtaOverlapCut(0.05),
386  fReco(0),
387  fFilterK0S(false),
388  fK0SLoss(0)
389 {
390  // DefineOutput(1,TList::Class());
391 }
392 //____________________________________________________________________
394  : AliAnalysisTaskSE(other),
395  fTracklets(0),
396  // fContainer(0),
397  fScaleDTheta(other.fScaleDTheta),
398  fMaxDelta(other.fMaxDelta),
399  fDPhiShift(other.fDPhiShift),
400  fDThetaWindow(other.fDThetaWindow),
401  fDPhiWindow(other.fDPhiWindow),
402  fPhiOverlapCut(other.fPhiOverlapCut),
403  fZEtaOverlapCut(other.fZEtaOverlapCut),
404  fReco(0),
405  fFilterK0S(other.fFilterK0S),
406  fK0SLoss(0)
407 {}
408 //____________________________________________________________________
410 {
411  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
412  if (!mgr) {
413  AliError("No analysis manager to connect to.");
414  return false;
415  }
416  AliAODHandler* ah =
417  dynamic_cast<AliAODHandler*>(mgr->GetOutputEventHandler());
418  if (!ah) {
419  AliError("No AOD output handler!");
420  return false;
421  }
422 
423  // Add to the manager
424  mgr->AddTask(this);
425 
426  // Always connect input
427  mgr->ConnectInput(this, 0, mgr->GetCommonInputContainer());
428 
429  return true;
430 }
431 //____________________________________________________________________
433 {
434  Printf("%s: %s", ClassName(), GetName());
435  Printf(" %22s: %d", "Scale by sin^2(theta)", fScaleDTheta);
436  Printf(" %22s: %f", "Delta phi shift", fDPhiShift);
437  Printf(" %22s: %f", "max Delta", fMaxDelta);
438  Printf(" %22s: %f", "Delta theta window", fDThetaWindow);
439  Printf(" %22s: %f", "Delta phi window", fDPhiWindow);
440  Printf(" %22s: %f", "phi overlap cut", fPhiOverlapCut);
441  Printf(" %22s: %f", "z-eta overlap cut", fZEtaOverlapCut);
442  Printf(" %22s: %s", "Filter K^0_S", fFilterK0S ? "yes" : "no");
443 }
444 
445 //____________________________________________________________________
447 {
448  return *this;
449 }
450 //____________________________________________________________________
452 {
453  // fContainer = new Container;
454  // fContainer->SetOwner();
455  // fContainer->SetName(GetName());
456  // PostData(1, fContainer);
457 
458  // Create container of tracklets
459  if (WorkerInit()) return;
460 
461  AliWarning("Failed to initialize task on worker, making zombie");
462  SetZombie();
463 
464 }
465 
466 //____________________________________________________________________
468 {
469  Printf("Now initialising CDB");
470  AliAnalysisManager* anaMgr = AliAnalysisManager::GetAnalysisManager();
471  if (!anaMgr) {
472  AliError("No manager defined!");
473  return false;
474  }
475  // Check if we have the CDB connect task, and if so, do nothing as
476  // we rely on that task set up things properly.
477  const char* cdbNames[] = {"CDBconnect", "cdb", 0 };
478  const char** ptr = cdbNames;
479  while (*ptr) {
480  AliAnalysisTask* cdbConnect = anaMgr->GetTask(*ptr);
481  if (cdbConnect && cdbConnect->IsA()->InheritsFrom("AliTaskCDBconnect")) {
482  AliInfoF("CDB-connect task (%s: %s) present, do nothing",
483  cdbConnect->ClassName(), *ptr);
484  return true;
485  }
486  ptr++;
487  }
488  // Otherwise, we need to do stuff ourselves
489  Printf("Get the CDB manager");
490  AliCDBManager* cdbMgr = AliCDBManager::Instance();
491  if (!cdbMgr) {
492  AliError("Failed to get instance of CDB manager");
493  return false;
494  }
495  Int_t refRun = GetCDBReferenceRun();
496  TString refUrl = GetCDBReferenceURL();
497  AliWarningF("Using reference CDB storage \"%s\" and run \"%d\"",
498  refUrl.Data(), refRun);
499  // Set a very particular default storage. Perhaps we can do this
500  // with specific storages instead! Depends on whether the
501  // reconstruction also uses CDB - probably does - in which case we
502  // need specific storages for that too.
503  cdbMgr->SetDefaultStorage(refUrl);
504  // Now load our geometry - from a LHC10h run
505  Printf("Get Geometry entry");
506  AliCDBEntry* cdbEnt = cdbMgr->Get("GRP/Geometry/Data", refRun);
507  if (!cdbEnt) {
508  AliErrorF("No geometry found from %d", refRun);
509  return false;
510  }
511  // Initialize the geometry manager
512  Printf("Set Geometry");
513  AliGeomManager::SetGeometry(static_cast<TGeoManager*>(cdbEnt->GetObject()));
514  // Now perform mis-alignment - again based on an LHC10h run!
515  Printf("Misalign geometry");
516  if (!AliGeomManager::ApplyAlignObjsToGeom("ITS",refRun,-1,-1)) {
517  AliErrorF("Failed to misalign geometry from %d", refRun);
518  return false;
519  }
520  return true;
521 }
522 //____________________________________________________________________
524 {
525  fTracklets = new TClonesArray(TrackletClassName());
526  // fTracklets->SetName(GetName());
527  TObject* obj = fTracklets;
528  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
529  AliAODHandler* ah =
530  dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
531  if (!ah) {
532  AliWarning("No AOD output handler set in analysis manager");
533  return false;
534  }
535  ah->AddBranch("TClonesArray", &obj);
536  // AliAODEvent* aod = ah->GetAOD();
537  // aod->Print();
538  // ah->GetTree()->Print();
539 
540  if (fFilterK0S) {
541  fK0SLoss = new TParameter<double>("k0Loss", 0);
542  ah->AddBranch("TParameter<double>", &fK0SLoss);
543  }
544  return true;
545 }
546 
547 //____________________________________________________________________
549 {
550  if (!InitCDB()) return false;
551  if (!InitBranch()) return false;
552 
553  // PostData(1, fContainer);
554 
555  return true;
556 }
557 //____________________________________________________________________
559 {
560  // Clear container of tracklets
561  if (!fTracklets) {
562  AliWarning("Tracklets container not initialized, init must have failed!");
563  return;
564  }
565  fTracklets->Clear();
566 
567  if (!ProcessEvent()) return;
568 
569  MarkForStore();
570 }
571 
572 //____________________________________________________________________
574 {
575  AliVEvent* event = FindEvent();
576  if (!event) return false;
577 
578 
579  TTree* clusters = FindClusters();
580  const AliVVertex* ip = FindIP(event);
581  Bool_t ret = false;
582  if (!ip) return false;
583  if (clusters) {
584  if (!HasGeometry()) ret = false;
585  if (ret && !HasField(event)) ret = false;
586  // If we have clusters, then try to reconstruct the event (again)
587  if (ret) ret = Reconstruct(clusters, ip);
588  }
589  else
590  // Otherwise, use the stored tracklets (max Delta = 1.5)
591  ret = ProcessTracklets(true, event->GetMultiplicity());
592 
593  CleanClusters(clusters);
594 
595  return ret;
596 }
597 
598 //____________________________________________________________________
600 {
601  if (!AliGeomManager::GetGeometry()) {
602  AliError("No geometry loaded, needed for reconstruction");
603  AliError("Add the AliTaskCDBconnect to the train");
604  return false;
605  }
606  return true;
607 }
608 //____________________________________________________________________
610 {
611  if (!TGeoGlobalMagField::Instance()->GetField() &&
612  !event->InitMagneticField()) {
613  AliWarning("Failed to initialize magnetic field");
614  return false;
615  }
616  return true;
617 }
618 
619 //____________________________________________________________________
621 {
622  AliVEvent* event = InputEvent();
623  if (!event) AliWarning("No event");
624  return event;
625 }
626 
627 //____________________________________________________________________
629 {
630  AliAnalysisManager* mgr = AliAnalysisManager::GetAnalysisManager();
631  AliVEventHandler* inh = mgr->GetInputEventHandler();
632  if (!inh->IsA()->InheritsFrom(AliESDInputHandlerRP::Class())) {
633  AliErrorF("Not the right kind of input handler: %s",
634  inh->ClassName());
635  return 0;
636  }
637  AliESDInputHandlerRP* rph = static_cast<AliESDInputHandlerRP*>(inh);
638  TTree* tree = rph->GetTreeR("ITS");
639  if (!tree) {
640  AliError("Tree of clusters (rec.points) not found");
641  return 0;
642  }
643  return FilterClusters(tree);
644 }
645 
646 //____________________________________________________________________
647 const AliVVertex* AliTrackletAODTask::FindIP(AliVEvent* event,
648  Double_t maxDispersion,
649  Double_t maxZError)
650 {
651  const AliVVertex* ip = event->GetPrimaryVertex();
652  if (ip->GetNContributors() <= 0) {
653  AliWarning("Not enough contributors for IP");
654  return 0;
655  }
656  // If this is from the Z vertexer, do some checks
657  if (ip->IsFromVertexerZ()) {
658  // Get covariance matrix
659  Double_t covar[6];
660  ip->GetCovarianceMatrix(covar);
661  Double_t sigmaZ = TMath::Sqrt(covar[5]);
662  if (sigmaZ >= maxZError) {
663  AliWarningF("IPz resolution = %f >= %f", sigmaZ, maxZError);
664  return 0;
665  }
666 
667  // If this IP doesn not derive from AliVertex, don't check dispersion.
668  if (ip->IsA()->InheritsFrom(AliVertex::Class())) {
669  const AliVertex* ipv = static_cast<const AliVertex*>(ip);
670  // Dispersion is the parameter used by the vertexer for finding the IP.
671  if (ipv->GetDispersion() >= maxDispersion) {
672  AliWarningF("IP dispersion = %f >= %f",
673  ipv->GetDispersion(), maxDispersion);
674  return 0;
675  }
676  }
677  }
678 
679 #if 0
680  // If we get here, we either have a full 3D vertex or track
681  // vertex, and we should check if it is in range
682  if (ip->GetZ() < fIPzAxis.GetXmin() || ip->GetZ() > fIPzAxis.GetXmax()) {
683  AliWarningF("IPz = %fcm out of range [%f,%f]cm",
684  ip->GetZ(), fIPzAxis.GetXmin(), fIPzAxis.GetXmax());
685  return 0;
686  }
687 #endif
688  // Good vertex, return it
689  return ip;
690 }
691 
692 //____________________________________________________________________
694 {
695 
696  AliAnalysisManager* am = AliAnalysisManager::GetAnalysisManager();
697  AliAODHandler* ah =
698  dynamic_cast<AliAODHandler*>(am->GetOutputEventHandler());
699  if (!ah) return false;
700  ah->SetFillAOD(kTRUE);
701  return true;
702 }
703 
704 //____________________________________________________________________
705 Bool_t AliTrackletAODTask::Reconstruct(TTree* clusters, const AliVVertex* ip)
706 {
707  // Make an IP array for the reconstructor
708  Float_t ipv[3];
709  ipv[0] = ip->GetX();
710  ipv[1] = ip->GetY();
711  ipv[2] = ip->GetZ();
712 
713  // Make the reconstructor
714  fReco = 0;
715  AliITSMultRecBg reco;
716  reco.SetCreateClustersCopy (true);
717  reco.SetScaleDThetaBySin2T (fScaleDTheta);
718  reco.SetNStdDev (fMaxDelta);
719  reco.SetPhiWindow (fDPhiWindow);
720  reco.SetThetaWindow (fDThetaWindow);
721  reco.SetPhiShift (fDPhiShift);
722  reco.SetRemoveClustersFromOverlaps(fPhiOverlapCut > 0 ||
723  fZEtaOverlapCut > 0);
724  reco.SetPhiOverlapCut (fPhiOverlapCut);
725  reco.SetZetaOverlapCut (fZEtaOverlapCut);
726  reco.SetHistOn (false);
727  reco.SetRecType (AliITSMultRecBg::kData);
728 
729  // Run normal reconstruction
730  reco.Run(clusters, ipv);
731 
732  // Save pointer to reco object - for MC combinatorics
733  fReco = &reco;
734  // And fill results into output branch
735  if (!ProcessTracklets(true, reco.GetMultiplicity())) return false;
736  fReco = 0;
737 
738  // Run again, but in injection mode
739  reco.SetRecType(AliITSMultRecBg::kBgInj);
740  reco.Run(clusters, ipv);
741 
742  // And fill results into output branch
743  if(!ProcessTracklets(false, reco.GetMultiplicity())) return false;
744 
745  return true;
746 }
747 
748 //____________________________________________________________________
750  AliVMultiplicity* vmult)
751 {
752  // Now loop over all tracklets. Since this method is only called
753  // once per "reconstruction" pass, we save a lot of computing time
754  // since we loop over the tracklets of each reconstruction exactly
755  // once.
756  // if (!mult->IsA()->InheritsFrom(AliMultiplicity::Class())) return false;
757  AliMultiplicity* mult = static_cast<AliMultiplicity*>(vmult);
758  Int_t nTracklets = mult->GetNumberOfTracklets();
759  for (Int_t trackletNo = 0; trackletNo < nTracklets; trackletNo++) {
760  if (!ProcessTracklet(normal,mult,trackletNo)) return false;
761  }
762  return true;
763 }
764 
765 //____________________________________________________________________
767 {
768  if (!fTracklets) return 0;
769  Int_t n = fTracklets->GetEntries();
770  return new((*fTracklets)[n]) AliAODTracklet;
771 }
772 //____________________________________________________________________
775  AliMultiplicity* mult,
776  Int_t no)
777 {
778  Double_t theta = mult->GetTheta (no);
779  Double_t phi = mult->GetPhi (no);
780  Double_t dTheta = mult->GetDeltaTheta(no);
781  Double_t dPhi = mult->GetDeltaPhi (no);
782  Double_t delta = mult->CalcDist (no);
783  AliAODTracklet* tracklet = MakeTracklet(normal);
784  tracklet->SetTheta (theta);
785  tracklet->SetPhi (phi);
786  tracklet->SetDTheta(dTheta);
787  tracklet->SetDPhi (dPhi);
788  tracklet->SetDelta (delta);
789  if (!normal) tracklet->SetInjection();
790 
791  return tracklet;
792 }
793 
794 /*********************************************************************
795  *
796  * Code for processing simulations
797  */
798 #ifndef __CINT__
799 #include <AliStack.h>
800 #include <AliMCEvent.h>
801 #include <AliGenEventHeader.h>
802 #include <TParticle.h>
803 #include <TParticlePDG.h>
804 #include <TDatabasePDG.h>
805 #include <AliITSgeomTGeo.h>
806 #include <TRandom.h>
807 #include <TROOT.h>
808 #else
809 class TParticlePDG;
810 class TParticle;
811 class AliStack;
812 #endif
813 //====================================================================
820 {
821 public:
831  AliTrackletAODMCTask(const char* name)
832  : AliTrackletAODTask(name)
833  {}
840  : AliTrackletAODTask(other)
841  {}
850  {
851  return *this;
852  }
853 private:
858  /*
859  * Initialize the worker
860  */
861  Bool_t WorkerInit();
867  virtual TTree* FilterClusters(TTree* t);
873  virtual void CleanClusters(TTree*& t);
880  virtual Bool_t ProcessEvent();
887  /* @} */
897  virtual const char* GetCDBReferenceURL() const
898  {
899  return "alien://Folder=/alice/simulation/2008/v4-15-Release/Residual";
900  }
901  /* @} */
911  virtual const char* TrackletClassName() const { return "AliAODMCTracklet"; }
919  virtual AliAODTracklet* MakeTracklet(Bool_t normal);
929  virtual AliAODTracklet* ProcessTracklet(Bool_t normal,
930  AliMultiplicity* mult,
931  Int_t no);
932  /* @} */
945  TParticle* FindPrimaryParent(Int_t label) const;
954  Int_t FindParent(Int_t label) const;
969  Int_t FindParents(Int_t label, TArrayI& fill, Int_t offset) const;
980  Int_t CommonParent(Int_t label, const TArrayI& fill) const;
981  /* @} */
982 
983 
985 };
986 
987 #define ADDP(P,M,L,Q) \
988  db->AddParticle("p" #P, "p" #P,M,false,L,Q,"Baryon",P); \
989  db->AddAntiParticle("p" #P "_bar", -P)
990 
991 //____________________________________________________________________
993 {
994  if (!AliTrackletAODTask::WorkerInit()) return false;
995 
996  // Register additional particles from EPOS-LHC
997  TDatabasePDG* db = TDatabasePDG::Instance();
998  db->GetParticle(2212);
999  db->AddParticle("deut","deut",1.876560,false,0.000000,3,"Baryon",1000010020);
1000  db->AddAntiParticle("deut_bar",-1000010020);
1001  db->AddParticle("trit","trit",2.816700,false,0.000000,3,"Baryon",1000010030);
1002  db->AddAntiParticle("trit_bar",-1000010030);
1003  db->AddParticle("alph","alph",3.755000,false,0.000000,6,"Baryon",1000020040);
1004  db->AddAntiParticle("alph_bar",-1000020040);
1005  ADDP(32224, 1.524000,0.145000,6);
1006  ADDP(12224, 1.816000,0.350000,6);
1007  ADDP(12222, 2.108000,0.350000,6);
1008  ADDP(12212, 1.525720,0.145000,3);
1009  ADDP(2124, 1.819440,0.150000,3);
1010  ADDP(32214, 2.113160,0.150000,3);
1011  ADDP(32212, 2.406880,0.145000,3);
1012  ADDP(12214, 2.700600,0.300000,3);
1013  ADDP(22124, 2.994320,0.150000,3);
1014  ADDP(12122, 3.288040,0.300000,3);
1015  ADDP(13222, 1.575200,0.080000,3);
1016  ADDP(23222, 1.768100,0.100000,3);
1017  ADDP(13226, 1.961000,0.170000,3);
1018  ADDP(12112, 1.524430,0.350000,0);
1019  ADDP(1214, 1.816860,0.150000,0);
1020  ADDP(32114, 2.109290,0.150000,0);
1021  ADDP(32112, 2.401720,0.145000,0);
1022  ADDP(12114, 2.694150,0.300000,0);
1023  ADDP(21214, 2.986579,0.150000,0);
1024  ADDP(11212, 3.279010,0.300000,0);
1025  ADDP(13122, 1.761000,0.040000,0);
1026  ADDP(3124, 1.950500,0.016000,0);
1027  ADDP(23122, 2.140000,0.090000,0);
1028  ADDP(13212, 2.329500,0.080000,0);
1029  ADDP(23212, 2.519000,0.100000,0);
1030  ADDP(43122, 2.708500,0.145000,0);
1031  ADDP(13216, 2.898000,0.170000,0);
1032  ADDP(31114, 1.524000,0.145000,-3);
1033  ADDP(11114, 1.816000,0.350000,-3);
1034  ADDP(11112, 2.108000,0.350000,-3);
1035  ADDP(13112, 1.577600,0.080000,-3);
1036  ADDP(23112, 1.767700,0.100000,-3);
1037  ADDP(13116, 1.957800,0.170000,-3);
1038  ADDP(9900110,0.000000,0.000000,0);
1039  ADDP(9900210,0.000000,0.000000,0);
1040  ADDP(9900220,0.000000,0.000000,0);
1041  ADDP(9900330,0.000000,0.000000,0);
1042  ADDP(9900440,0.000000,0.000000,0);
1043  ADDP(9902210,0.000000,0.000000,0);
1044  ADDP(9902110,0.000000,0.000000,0);
1045  ADDP(88, 0.000000,0.000000,0);
1046  ADDP(90, 0.000000,0.000000,0);
1047  ADDP(990, 0.000000,0.000000,0);
1048  ADDP(99999, 0.000000,0.000000,0);
1049 
1050  return true;
1051 }
1052 
1053 //____________________________________________________________________
1055 {
1056  if (fK0SLoss) fK0SLoss->SetVal(0);
1057  if (!t || !fFilterK0S) return t;
1058 
1059  const Double_t k0s = 1.52233299626516083e+00; // 310 - K^0_S weight
1060  const Double_t kpm = (1.43744204476109627e+00*
1061  9.82150320171071400e-01); // 321 - K^{+/-}
1062  const Double_t lam = 2.75002089647900005e+00; // 3122 - lambda
1063  const Double_t sig = 2.75002089647899961e+00; // 3212 - sigma
1064  const Double_t xi = 3.24109605656453548e+00; // 3322 - Xi
1065  TDirectory* savDir = gDirectory;
1066  gDirectory = gROOT;
1067  TTree* copy = new TTree("TreeR", "TreeR");
1068  TClonesArray* in = new TClonesArray("AliITSRecPoint");
1069  TClonesArray* out = new TClonesArray("AliITSRecPoint");
1070  Int_t outN = 0;
1071  copy->Branch("ITSRecPoints", &out);
1072  copy->SetAutoFlush(0);// Keep in memory
1073  t->SetBranchAddress("ITSRecPoints", &in);
1074 
1075  // Printf("Filtering clusters from strange particles");
1076 
1077  Int_t min1 = AliITSgeomTGeo::GetModuleIndex(1,1,1);
1078  Int_t max1 = AliITSgeomTGeo::GetModuleIndex(2,1,1);
1079  Int_t min2 = AliITSgeomTGeo::GetModuleIndex(2,1,1);
1080  Int_t max2 = AliITSgeomTGeo::GetModuleIndex(3,1,1);
1081  Int_t nTotal = 0;
1082  Int_t nKept = 0;
1083 
1084  // Loop over the modules of the SPD
1085  for (Int_t i = 0; i < max2; i++) {
1086  in->Clear();
1087  out->Clear();
1088  outN = 0;
1089 
1090  // Read in module data
1091  t->GetEntry(i);
1092 
1093  // Loop over all clusters in the module
1094  Int_t inN = in->GetEntries();
1095  for (Int_t j = 0; j < inN; j++) {
1096  AliITSRecPoint* inCl = static_cast<AliITSRecPoint*>(in->At(j));
1097  if (!inCl) continue;
1098 
1099  // Loop over labels in the module
1100  Double_t weight = 1;
1101  for (Int_t k = 0; k < 3; k++) {
1102  Int_t label = inCl->GetLabel(k);
1103  if (label <= 0) continue;
1104 
1105  // Check primary parent particle type
1106  TParticle* parent = FindPrimaryParent(label);
1107  if (!parent) continue;
1108 
1109  switch (TMath::Abs(parent->GetPdgCode())) {
1110  case 310: weight *= k0s; break;
1111  case 321: weight *= kpm; break;
1112  case 3122: weight *= lam; break;
1113  case 3212: weight *= sig; break;
1114  case 3322: weight *= xi; break;
1115  }
1116  if (weight > 1)
1117  Printf("Cluster from %+5d: %7.5f", parent->GetPdgCode(), weight);
1118  // Printf("Parent %d is a K^0_S: %d", k+1, parent->GetPdgCode());
1119  // Randomly remove clusters from the right kind of primary
1120  // parent particles
1121  }
1122  if (weight > 1) {
1123  Double_t chance = 1 - 1 / weight; // chance is 1 minus inverse
1124  nTotal++;
1125  if (gRandom->Uniform() < chance) continue;
1126  nKept++;
1127  }
1128  new ((*out)[outN++]) AliITSRecPoint(*inCl);
1129  }
1130  // Printf("Kept %d out of %d clusters", outN, inN);
1131  copy->Fill();
1132  }
1133  if (nTotal > 0) {
1134  Double_t loss = 1-Float_t(nKept)/nTotal;
1135  Printf("Kept %d out of %d clusters from strange primaries (%4.1f%%)",
1136  nKept, nTotal, 100.*loss);
1137  if (fK0SLoss) fK0SLoss->SetVal(loss);
1138  }
1139  savDir->cd();
1140  return copy;
1141 }
1142 
1143 //____________________________________________________________________
1145 {
1146  if (!t || !fFilterK0S) return;
1147 
1148  delete t;
1149  t = 0;
1150 }
1151 
1152 //____________________________________________________________________
1154 {
1155  if (!AliTrackletAODTask::ProcessEvent()) return false;
1156 
1157  // create generated "tracklets"
1158  if (!ProcessGenerated()) return false;
1159 
1160  return true;
1161 }
1162 
1163 //____________________________________________________________________
1165 {
1166  AliStack* stack = MCEvent()->Stack();
1167  Int_t nTracks = stack->GetNtrack();
1168  for (Int_t trackNo = nTracks; trackNo--; ) {
1169  if (!stack->IsPhysicalPrimary(trackNo)) {
1170  // AliWarningF("Track # %6d not a primary", trackNo);
1171  // Not a primary, go on
1172  continue;
1173  }
1174  // Get the particle
1175  TParticle* particle = stack->Particle(trackNo);
1176  if (!particle) {
1177  AliWarningF("No particle found for track # %d", trackNo);
1178  continue;
1179  }
1180 
1181  TParticlePDG* pdg = particle->GetPDG();
1182  if (!pdg) {
1183  ::Warning("GetPDG", "Unknown PDG code: %d", particle->GetPdgCode());
1184  // continue; // Unknown particle
1185  }
1186  // if (pdg->Charge() == 0) {
1187  // Uncharged - don't care
1188  // continue;
1189  // }
1190 
1191  // Get theta
1192  Double_t theta = particle->Theta();
1193  // Check for beam-like particle
1194  if (theta < 1e-6 || TMath::Abs(theta-TMath::Pi()) < 1e-6) {
1195  AliWarningF("Track # %6d is beam-like (%f)", trackNo,
1196  TMath::RadToDeg()*theta);
1197  continue;
1198  }
1199  Double_t eta = -TMath::Log(TMath::ATan(theta/2));
1200  // If the pseudorapidity is way beyond the SPD acceptance, do not
1201  // write down this generated "tracklet" - to save space.
1202  if (TMath::Abs(eta) > 3) continue;
1203 
1204  Double_t phi = particle->Phi();
1205  AliAODTracklet* tracklet = MakeTracklet(false);
1206  AliAODMCTracklet* mc = static_cast<AliAODMCTracklet*>(tracklet);
1207  mc->SetGenerated();
1208  mc->SetTheta(theta);
1209  mc->SetPhi(phi);
1210  mc->SetDTheta(0);
1211  mc->SetDPhi (0);
1212  mc->SetDelta (0);
1213  mc->SetParentPdg(particle->GetPdgCode());
1214  mc->SetParentPt (particle->Pt());
1215  if (pdg && pdg->Charge() == 0) mc->SetNeutral();
1216  }
1217  return true;
1218 }
1219 
1220 //____________________________________________________________________
1222 {
1223  if (!fTracklets) return 0;
1224  Int_t n = fTracklets->GetEntries();
1225  return new((*fTracklets)[n]) AliAODMCTracklet;
1226 }
1227 //____________________________________________________________________
1230  AliMultiplicity* mult,
1231  Int_t no)
1232 {
1233  AliAODTracklet* tracklet =
1234  AliTrackletAODTask::ProcessTracklet(normal,mult,no);
1235  if (!normal) return tracklet;
1236 
1237  AliAODMCTracklet* mc = static_cast<AliAODMCTracklet*>(tracklet);
1238  Int_t label0 = mult->GetLabel(no, 0);
1239  Int_t label1 = mult->GetLabel(no, 1);
1240  TParticle* parent0 = FindPrimaryParent(label0);
1241  if (parent0) {
1242  mc->SetParentPdg(parent0->GetPdgCode());
1243  mc->SetParentPt (parent0->Pt());
1244  }
1245  if (label0 != label1) {
1246  TParticle* parent1 = FindPrimaryParent(label1);
1247  if (parent1) {
1248  mc->SetParentPdg(parent1->GetPdgCode(), true);
1249  mc->SetParentPt (parent1->Pt(), true);
1250  }
1251  mc->SetCombinatorics();
1252  // Here, we could track back in the cluster labels to see if we
1253  // have the same ultimate mother of both clusters. Since this
1254  // isn't really used, we do not do that
1255  if (fReco) {
1256  TArrayI parents(50); // At most 50 levels deep
1257  // Tracklet parameters
1258  Float_t* ftrack = fReco->GetTracklet(no);
1259  // Cluster identifiers
1260  Int_t clus0 = Int_t(ftrack[AliITSMultReconstructor::kClID1]);
1261  Int_t clus1 = Int_t(ftrack[AliITSMultReconstructor::kClID2]);
1262  // Cluster labelsx
1263  Float_t* fclus0 = (fReco->GetClusterOfLayer(0,clus0) +
1264  AliITSMultReconstructor::kClMC0);
1265  Float_t* fclus1 = (fReco->GetClusterOfLayer(1,clus1) +
1266  AliITSMultReconstructor::kClMC0);
1267  // Loop over three inner layer cluster labels
1268  Int_t offset = 0;
1269  for (Int_t i = 0; i < 3; i++)
1270  offset = FindParents(Int_t(fclus0[i]), parents, offset);
1271  // Loop over three outer layer cluster labels
1272  Bool_t distinct = true;
1273  for (Int_t i = 0; i < 3; i++) {
1274  Float_t flbl = fclus1[i];
1275  // Be careful not to get overflows on conversions
1276  if (flbl > nextafter(INT_MAX, 0) || flbl < nextafter(INT_MIN, 0))
1277  continue; // Would overflow
1278  if (CommonParent(Int_t(flbl), parents)) {
1279  distinct = false;
1280  // We break out as soon as we find a common parent.
1281  break;
1282  }
1283  } // loop over outer layer
1284  if (distinct) mc->SetDistinct();
1285  } // if (fReco)
1286  }
1287  else {
1288  if (!MCEvent()->Stack()->IsPhysicalPrimary(label0))
1289  mc->SetSecondary();
1290  }
1291  return tracklet;
1292 }
1293 
1294 //____________________________________________________________________
1296 {
1297  AliStack* stack = MCEvent()->Stack();
1298  Int_t nTracks = stack->GetNtrack();
1299  TParticle* particle = stack->Particle(label);
1300  if (!particle) return -1;
1301  Int_t ret = particle->GetFirstMother();
1302  if (ret > nTracks) return -1;
1303  return ret;
1304 }
1305 
1306 //____________________________________________________________________
1308  TArrayI& fill,
1309  Int_t offset) const
1310 {
1311  // If offset is negative, then we don't really store the labels, but
1312  // only check against those already stored.
1313  Int_t i = offset;
1314  Int_t lbl = label;
1315  while (i < fill.GetSize()-1 && lbl >= 1) {
1316  fill[i] = lbl;
1317  i++;
1318  lbl = FindParent(lbl);
1319  }
1320  // If we get here, and we're checking, that means we did not find a
1321  // common ancestor, so we return 0 (false). If we're not checking,
1322  // then we return the number of elements assigned in the passed
1323  // cache.
1324  return i;
1325 }
1326 
1327 //____________________________________________________________________
1329  const TArrayI& fill) const
1330 {
1331  Int_t i = 0;
1332  Int_t lbl = label;
1333  while (i < fill.GetSize()-1 && lbl >= 1) {
1334  // If we're checking, just see if we have the label in the list
1335  // already. If we do, then return the index
1336  for (Int_t j = 0; j < fill.GetSize(); j++) {
1337  if (fill[j] == lbl) return j;
1338  }
1339  i++;
1340  lbl = FindParent(lbl);
1341  }
1342  // If we get here, and we're checking, that means we did not find a
1343  // common ancestor, so we return 0 (false). If we're not checking,
1344  // then we return the number of elements assigned in the passed
1345  // cache.
1346  return false;
1347 }
1348 
1349 
1350 //____________________________________________________________________
1352 {
1353  AliStack* stack = MCEvent()->Stack();
1354  Int_t nTracks = stack->GetNtrack();
1355  TParticle* particle = stack->Particle(label);
1356  Int_t trackNo = label;
1357 #if 1
1358  while (!stack->IsPhysicalPrimary(trackNo)) {
1359  trackNo = particle->GetFirstMother();
1360  // If we have hit the top
1361  if (trackNo < 0) return 0;
1362  // Partice first next iteration
1363  particle = stack->Particle(trackNo);
1364  }
1365 
1366 #else
1367  Int_t parentID = particle->GetFirstMother();
1368  while (parentID >= 0 && parentID < nTracks) {
1369  particle = stack->Particle(parentID);
1370  parentID = particle->GetFirstMother();
1371  }
1372 #endif
1373  return particle;
1374 }
1375 //====================================================================
1377 {
1378  AliAnalysisManager *mgr = AliAnalysisManager::GetAnalysisManager();
1379  if (!mgr) {
1380  ::Error("Create","No analysis manager to connect to.");
1381  return 0;
1382  }
1383  AliAODHandler* ah =
1384  dynamic_cast<AliAODHandler*>(mgr->GetOutputEventHandler());
1385  if (!ah) {
1386  ::Error("Create","No AOD output handler!");
1387  return 0;
1388  }
1389  Bool_t mc = mgr->GetMCtruthEventHandler() != 0;
1390  AliTrackletAODTask* ret = 0;
1391  if (mc) ret = new AliTrackletAODMCTask("MidRapidityMC");
1392  else ret = new AliTrackletAODTask("MidRapidity");
1393  if (ret) ret->Connect();
1394 
1395  return ret;
1396 }
1397 
1398 
1399 //====================================================================
1400 
1401 //
1402 // EOF
1403 //
1404 
1405 
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 *)
TParameter< double > * fK0SLoss
Real_t GetTheta() const
virtual const char * GetCDBReferenceURL() 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
virtual void SetFilterK0S(Bool_t filter)
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)
virtual Bool_t ProcessEvent()
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 const char * GetCDBReferenceURL() const
AliTrackletAODMCTask(const AliTrackletAODMCTask &other)
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)
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 *)