Corrected protection.
[u/mrichter/AliRoot.git] / STEER / AliAODEvent.cxx
index d0d1497..65ed50e 100644 (file)
@@ -30,6 +30,7 @@
 #include "AliAODEvent.h"
 #include "AliAODHeader.h"
 #include "AliAODTrack.h"
+#include "AliAODDimuon.h"
 
 ClassImp(AliAODEvent)
 
@@ -45,7 +46,8 @@ ClassImp(AliAODEvent)
                                                      "phosCells",
                                                      "caloClusters",
                                                      "fmdClusters",
-                                                     "pmdClusters"
+                                                     "pmdClusters",
+                                                     "dimuons"
                                                      
 };
 //______________________________________________________________________________
@@ -65,7 +67,8 @@ AliAODEvent::AliAODEvent() :
   fPhosCells(0),
   fCaloClusters(0),
   fFmdClusters(0),
-  fPmdClusters(0)
+  fPmdClusters(0),
+  fDimuons(0)
 {
   // default constructor
 }
@@ -74,7 +77,7 @@ AliAODEvent::AliAODEvent() :
 AliAODEvent::AliAODEvent(const AliAODEvent& aod):
   AliVEvent(aod),
   fAODObjects(new TList()),
-  fAODFolder(new TFolder()),
+  fAODFolder(0),
   fConnected(kFALSE),
   fHeader(new AliAODHeader(*aod.fHeader)),
   fTracks(new TClonesArray(*aod.fTracks)),
@@ -87,7 +90,8 @@ AliAODEvent::AliAODEvent(const AliAODEvent& aod):
   fPhosCells(new AliAODCaloCells(*aod.fPhosCells)),
   fCaloClusters(new TClonesArray(*aod.fCaloClusters)),
   fFmdClusters(new TClonesArray(*aod.fFmdClusters)),
-  fPmdClusters(new TClonesArray(*aod.fPmdClusters))
+  fPmdClusters(new TClonesArray(*aod.fPmdClusters)),
+  fDimuons(new TClonesArray(*aod.fDimuons))
 {
   // Copy constructor
   AddObject(fHeader);
@@ -102,6 +106,7 @@ AliAODEvent::AliAODEvent(const AliAODEvent& aod):
   AddObject(fCaloClusters);
   AddObject(fFmdClusters);
   AddObject(fPmdClusters);
+  AddObject(fDimuons);
   fConnected = aod.fConnected;
   GetStdContent();
   CreateStdFolders();
@@ -214,6 +219,7 @@ AliAODEvent::~AliAODEvent()
     }
 
     delete fAODFolder;
+    fAODFolder = 0;
 }
 
 //______________________________________________________________________________
@@ -235,7 +241,7 @@ void AliAODEvent::RemoveObject(TObject* obj)
 }
 
 //______________________________________________________________________________
-TObject *AliAODEvent::FindListObject(const char *objName)
+TObject *AliAODEvent::FindListObject(const char *objName) const
 {
   // Return the pointer to the object with the given name.
 
@@ -260,6 +266,7 @@ void AliAODEvent::CreateStdContent()
   AddObject(new TClonesArray("AliAODCaloCluster", 0));
   AddObject(new TClonesArray("AliAODFmdCluster", 0));
   AddObject(new TClonesArray("AliAODPmdCluster", 0));
+  AddObject(new TClonesArray("AliAODDimuon", 0));
   // set names
   SetStdNames();
 
@@ -308,6 +315,7 @@ void AliAODEvent::SetStdNames()
 void AliAODEvent::CreateStdFolders()
 {
     // Create the standard folder structure
+  if(fAODFolder)delete fAODFolder;
     fAODFolder = gROOT->GetRootFolder()->AddFolder("AOD", "AOD");
     if(fAODObjects->GetEntries()==kAODListN){
        for(int i = 0;i < fAODObjects->GetEntries();i++){
@@ -341,6 +349,7 @@ void AliAODEvent::GetStdContent()
   fCaloClusters  = (TClonesArray*)fAODObjects->FindObject("caloClusters");
   fFmdClusters   = (TClonesArray*)fAODObjects->FindObject("fmdClusters");
   fPmdClusters   = (TClonesArray*)fAODObjects->FindObject("pmdClusters");
+  fDimuons       = (TClonesArray*)fAODObjects->FindObject("dimuons");
 }
 
 //______________________________________________________________________________
@@ -351,7 +360,8 @@ void AliAODEvent::ResetStd(Int_t trkArrSize,
                           Int_t jetSize, 
                           Int_t caloClusSize, 
                           Int_t fmdClusSize, 
-                          Int_t pmdClusSize
+                          Int_t pmdClusSize,
+                          Int_t dimuonArrSize
                           )
 {
   // deletes content of standard arrays and resets size 
@@ -363,7 +373,7 @@ void AliAODEvent::ResetStd(Int_t trkArrSize,
   fVertices->Delete();
   if (vtxArrSize > fVertices->GetSize()) 
     fVertices->Expand(vtxArrSize);
-        
+
   fV0s->Delete();
   if (v0ArrSize > fV0s->GetSize()) 
     fV0s->Expand(v0ArrSize);
@@ -387,6 +397,10 @@ void AliAODEvent::ResetStd(Int_t trkArrSize,
   fPmdClusters->Delete();
   if (pmdClusSize > fPmdClusters->GetSize()) 
     fPmdClusters->Expand(pmdClusSize);
+    
+  fDimuons->Delete();
+  if (dimuonArrSize > fDimuons->GetSize()) 
+    fDimuons->Expand(dimuonArrSize);
 
   // Reset the tracklets
   fTracklets->DeleteContainer();
@@ -410,6 +424,7 @@ void AliAODEvent::ClearStd()
   fCaloClusters  ->Delete();
   fFmdClusters   ->Clear();
   fPmdClusters   ->Clear();
+  fDimuons       ->Clear();
 }
 
 //_________________________________________________________________
@@ -421,9 +436,9 @@ Int_t AliAODEvent::GetPHOSClusters(TRefArray *clusters) const
   
   AliAODCaloCluster *cl = 0;
   Bool_t first = kTRUE;
-  for (Int_t i = 0; i < GetNCaloClusters() ; i++) {
+  for (Int_t i = 0; i < GetNumberOfCaloClusters() ; i++) {
     if ( (cl = GetCaloCluster(i)) ) {
-      if (cl->IsPHOSCluster()){
+      if (cl->IsPHOS()){
        if(first) {
          new (clusters) TRefArray(TProcessID::GetProcessWithUID(cl)); 
          first=kFALSE;
@@ -444,9 +459,9 @@ Int_t AliAODEvent::GetEMCALClusters(TRefArray *clusters) const
   clusters->Clear();
   AliAODCaloCluster *cl = 0;
   Bool_t first = kTRUE;
-  for (Int_t i = 0; i < GetNCaloClusters(); i++) {
+  for (Int_t i = 0; i < GetNumberOfCaloClusters(); i++) {
     if ( (cl = GetCaloCluster(i)) ) {
-      if (cl->IsEMCALCluster()){
+      if (cl->IsEMCAL()){
        if(first) {
          new (clusters) TRefArray(TProcessID::GetProcessWithUID(cl)); 
          first=kFALSE;
@@ -469,7 +484,8 @@ Int_t AliAODEvent::GetMuonTracks(TRefArray *muonTracks) const
 
   AliAODTrack *track = 0;
   for (Int_t iTrack = 0; iTrack < GetNTracks(); iTrack++) {
-    if ((track = GetTrack(iTrack))->IsMuonTrack()) {
+    track = GetTrack(iTrack);
+    if (track->IsMuonTrack()) {
       muonTracks->Add(track);
     }
   }
@@ -478,103 +494,118 @@ Int_t AliAODEvent::GetMuonTracks(TRefArray *muonTracks) const
 }
 
 
+//______________________________________________________________________________
+Int_t AliAODEvent::GetNumberOfMuonTracks() const
+{
+  // get number of muon tracks
+  Int_t nMuonTracks=0;
+  AliAODTrack *track = 0;
+  for (Int_t iTrack = 0; iTrack < GetNTracks(); iTrack++) {
+    if ((track = GetTrack(iTrack))->IsMuonTrack()) {
+       nMuonTracks++;
+    }
+  }
+  
+  return nMuonTracks;
+}
+
 void AliAODEvent::ReadFromTree(TTree *tree, Option_t* opt /*= ""*/)
 {
-  // Connects aod event to tree
+    // Connects aod event to tree
   
   if(!tree){
     Printf("%s %d AliAODEvent::ReadFromTree() Zero Pointer to Tree \n",(char*)__FILE__,__LINE__);
     return;
   }
-  // load the TTree
+    // load the TTree
   if(!tree->GetTree())tree->LoadTree(0);
-
-  // Try to find AliAODEvent
+  
+    // Try to find AliAODEvent
   AliAODEvent *aodEvent = 0;
   aodEvent = (AliAODEvent*)tree->GetTree()->GetUserInfo()->FindObject("AliAODEvent");
   if(aodEvent){
-    // Check if already connected to tree
+      // Check if already connected to tree
     TList* connectedList = (TList*) (tree->GetUserInfo()->FindObject("AODObjectsConnectedToTree"));
     if (connectedList && (strcmp(opt, "reconnect"))) {
-       // If connected use the connected list of objects
+        // If connected use the connected list of objects
        fAODObjects->Delete();
        fAODObjects = connectedList;
        GetStdContent(); 
        fConnected = kTRUE;
        return;
     } 
-    // Connect to tree
-    // prevent a memory leak when reading back the TList
-    // if (!(strcmp(opt, "reconnect"))) fAODObjects->Delete();
-
-    // create a new TList from the UserInfo TList... 
-    // copy constructor does not work...
+      // Connect to tree
+      // prevent a memory leak when reading back the TList
+//      if (!(strcmp(opt, "reconnect"))) fAODObjects->Delete();
+    
+      // create a new TList from the UserInfo TList... 
+      // copy constructor does not work...
     fAODObjects = (TList*)(aodEvent->GetList()->Clone());
     fAODObjects->SetOwner(kTRUE);
     if(fAODObjects->GetEntries()<kAODListN){
       printf("%s %d AliAODEvent::ReadFromTree() TList contains less than the standard contents %d < %d \n",
-            (char*)__FILE__,__LINE__,fAODObjects->GetEntries(),kAODListN);
+             (char*)__FILE__,__LINE__,fAODObjects->GetEntries(),kAODListN);
     }
-    //
-    // Let's find out whether we have friends
+      //
+      // Let's find out whether we have friends
     TList* friendL = tree->GetTree()->GetListOfFriends();
     if (friendL) 
     {
-       TIter next(friendL);
-       TFriendElement* fe;
-       while ((fe = (TFriendElement*)next())){
-           aodEvent = (AliAODEvent*)(fe->GetTree()->GetUserInfo()->FindObject("AliAODEvent"));
-           if (!aodEvent) {
-               printf("No UserInfo on tree \n");
-           } else {
-
-               TList* objL = (TList*)(aodEvent->GetList()->Clone());
-               printf("Get list of object from tree %d !!\n", objL->GetEntries());
-               TIter nextobject(objL);
-               TObject* obj =  0;
-               while((obj = nextobject()))
-               {
-                   printf("Adding object from friend %s !\n", obj->GetName());
-                   fAODObjects->Add(obj);
-               } // object "branch" loop
-           } // has userinfo  
-       } // friend loop
+      TIter next(friendL);
+      TFriendElement* fe;
+      while ((fe = (TFriendElement*)next())){
+        aodEvent = (AliAODEvent*)(fe->GetTree()->GetUserInfo()->FindObject("AliAODEvent"));
+        if (!aodEvent) {
+          printf("No UserInfo on tree \n");
+        } else {
+          
+          TList* objL = (TList*)(aodEvent->GetList()->Clone());
+          printf("Get list of object from tree %d !!\n", objL->GetEntries());
+          TIter nextobject(objL);
+          TObject* obj =  0;
+          while((obj = nextobject()))
+          {
+            printf("Adding object from friend %s !\n", obj->GetName());
+            fAODObjects->Add(obj);
+          } // object "branch" loop
+        } // has userinfo  
+      } // friend loop
     } // has friends   
-           
-
-// set the branch addresses
+    
+    
+      // set the branch addresses
     TIter next(fAODObjects);
     TNamed *el;
     while((el=(TNamed*)next())){
       TString bname(el->GetName());
-      // check if branch exists under this Name
+        // check if branch exists under this Name
       TBranch *br = tree->GetTree()->GetBranch(bname.Data());
       if(br){
-       tree->SetBranchAddress(bname.Data(),fAODObjects->GetObjectRef(el));
+        tree->SetBranchAddress(bname.Data(),fAODObjects->GetObjectRef(el));
       } else {
-         br = tree->GetBranch(Form("%s.",bname.Data()));
-         if(br){
-             tree->SetBranchAddress(Form("%s.",bname.Data()),fAODObjects->GetObjectRef(el));
-         }
-         else{
-             printf("%s %d AliAODEvent::ReadFromTree() No Branch found with Name %s. \n",
-                    (char*)__FILE__,__LINE__,bname.Data());
-         }     
+        br = tree->GetBranch(Form("%s.",bname.Data()));
+        if(br){
+          tree->SetBranchAddress(Form("%s.",bname.Data()),fAODObjects->GetObjectRef(el));
+        }
+        else{
+          printf("%s %d AliAODEvent::ReadFromTree() No Branch found with Name %s. \n",
+                 (char*)__FILE__,__LINE__,bname.Data());
+        }      
       }
     }
     GetStdContent();
-    // when reading back we are not owner of the list 
-    // must not delete it
+      // when reading back we are not owner of the list 
+      // must not delete it
     fAODObjects->SetOwner(kTRUE);
     fAODObjects->SetName("AODObjectsConnectedToTree");
-    // we are not owner of the list objects 
-    // must not delete it
+      // we are not owner of the list objects 
+      // must not delete it
     tree->GetUserInfo()->Add(fAODObjects);
     fConnected = kTRUE;
   }// no aodEvent
   else {
-    // we can't get the list from the user data, create standard content
-    // and set it by hand
+      // we can't get the list from the user data, create standard content
+      // and set it by hand
     CreateStdContent();
     TIter next(fAODObjects);
     TNamed *el;
@@ -583,11 +614,123 @@ void AliAODEvent::ReadFromTree(TTree *tree, Option_t* opt /*= ""*/)
       tree->SetBranchAddress(bname.Data(),fAODObjects->GetObjectRef(el));
     }
     GetStdContent();
-    // when reading back we are not owner of the list 
-    // must not delete it
+      // when reading back we are not owner of the list 
+      // must not delete it
     fAODObjects->SetOwner(kTRUE);
   }
 }
+  //______________________________________________________________________________
+Int_t  AliAODEvent::GetNumberOfPileupVerticesSPD() const{
+  // count number of SPD pileup vertices
+  Int_t nVertices=GetNumberOfVertices();
+  Int_t nPileupVertices=0;
+  for(Int_t iVert=0; iVert<nVertices; iVert++){
+    AliAODVertex *v=GetVertex(iVert);
+    if(v->GetType()==AliAODVertex::kPileupSPD) nPileupVertices++;
+  }
+  return nPileupVertices;
+}
+//______________________________________________________________________________
+Int_t  AliAODEvent::GetNumberOfPileupVerticesTracks() const{
+  // count number of track pileup vertices
+  Int_t nVertices=GetNumberOfVertices();
+  Int_t nPileupVertices=0;
+  for(Int_t iVert=0; iVert<nVertices; iVert++){
+    AliAODVertex *v=GetVertex(iVert);
+    if(v->GetType()==AliAODVertex::kPileupTracks) nPileupVertices++;
+  }
+  return nPileupVertices;
+}
+//______________________________________________________________________________
+AliAODVertex* AliAODEvent::GetPrimaryVertexSPD() const{
+  //
+  Int_t nVertices=GetNumberOfVertices();
+  for(Int_t iVert=0; iVert<nVertices; iVert++){
+    AliAODVertex *v=GetVertex(iVert);
+    if(v->GetType()==AliAODVertex::kMainSPD) return v;
+  }
+  return 0;
+}
+//______________________________________________________________________________
+AliAODVertex* AliAODEvent::GetPileupVertexSPD(Int_t iV) const{
+  //
+  Int_t nVertices=GetNumberOfVertices();
+  Int_t counter=0;
+  for(Int_t iVert=0; iVert<nVertices; iVert++){
+    AliAODVertex *v=GetVertex(iVert);
+    if(v->GetType()==AliAODVertex::kPileupSPD){
+      if(counter==iV) return v;
+      ++counter;
+    }
+  }
+  return 0;
+}
+//______________________________________________________________________________
+AliAODVertex* AliAODEvent::GetPileupVertexTracks(Int_t iV) const{
+  //
+  Int_t nVertices=GetNumberOfVertices();
+  Int_t counter=0;
+  for(Int_t iVert=0; iVert<nVertices; iVert++){
+    AliAODVertex *v=GetVertex(iVert);
+    if(v->GetType()==AliAODVertex::kPileupTracks){
+      if(counter==iV) return v;
+      ++counter;
+    }
+  }
+  return 0;
+}
+//______________________________________________________________________________
+Bool_t  AliAODEvent::IsPileupFromSPD(Int_t minContributors, 
+                                    Double_t minZdist, 
+                                    Double_t nSigmaZdist, 
+                                    Double_t nSigmaDiamXY, 
+                                    Double_t nSigmaDiamZ) const{
+  //
+  // This function checks if there was a pile up
+  // reconstructed with SPD
+  //
+  AliAODVertex *mainV=GetPrimaryVertexSPD();
+  if(!mainV) return kFALSE;
+  Int_t nc1=mainV->GetNContributors();
+  if(nc1<1) return kFALSE;
+  Int_t nPileVert=GetNumberOfPileupVerticesSPD();
+  if(nPileVert==0) return kFALSE;
+  Int_t nVertices=GetNumberOfVertices();
+  
+  for(Int_t iVert=0; iVert<nVertices; iVert++){
+    AliAODVertex *pv=GetVertex(iVert);
+    if(pv->GetType()!=AliAODVertex::kPileupSPD) continue;
+    Int_t nc2=pv->GetNContributors();
+    if(nc2>=minContributors){
+      Double_t z1=mainV->GetZ();
+      Double_t z2=pv->GetZ();
+      Double_t distZ=TMath::Abs(z2-z1);
+      Double_t distZdiam=TMath::Abs(z2-GetDiamondZ());
+      Double_t cutZdiam=nSigmaDiamZ*TMath::Sqrt(GetSigma2DiamondZ());
+      if(GetSigma2DiamondZ()<0.0001)cutZdiam=99999.; //protection for missing z diamond information
+      if(distZ>minZdist && distZdiam<cutZdiam){
+       Double_t x2=pv->GetX();
+       Double_t y2=pv->GetY();
+       Double_t distXdiam=TMath::Abs(x2-GetDiamondX());
+       Double_t distYdiam=TMath::Abs(y2-GetDiamondY());
+       Double_t cov1[6],cov2[6];       
+       mainV->GetCovarianceMatrix(cov1);
+       pv->GetCovarianceMatrix(cov2);
+       Double_t errxDist=TMath::Sqrt(cov2[0]+GetSigma2DiamondX());
+       Double_t erryDist=TMath::Sqrt(cov2[2]+GetSigma2DiamondY());
+       Double_t errzDist=TMath::Sqrt(cov1[5]+cov2[5]);
+       Double_t cutXdiam=nSigmaDiamXY*errxDist;
+       if(GetSigma2DiamondX()<0.0001)cutXdiam=99999.; //protection for missing diamond information
+       Double_t cutYdiam=nSigmaDiamXY*erryDist;
+       if(GetSigma2DiamondY()<0.0001)cutYdiam=99999.; //protection for missing diamond information
+       if( (distXdiam<cutXdiam) && (distYdiam<cutYdiam) && (distZ>nSigmaZdist*errzDist) ){
+         return kTRUE;
+       }
+      }
+    }
+  }
+  return kFALSE;
+}
 
 //______________________________________________________________________________
 void AliAODEvent::Print(Option_t *) const
@@ -615,3 +758,11 @@ void AliAODEvent::AssignIDtoCollection(TCollection* col)
     while ((obj = next()))
        TProcessID::AssignID(obj);
 }
+
+Bool_t AliAODEvent::IsPileupFromSPDInMultBins() const {
+    Int_t nTracklets=GetTracklets()->GetNumberOfTracklets();
+    if(nTracklets<20) return IsPileupFromSPD(3,0.8);
+    else if(nTracklets<50) return IsPileupFromSPD(4,0.8);
+    else return IsPileupFromSPD(5,0.8);
+}
+