[u/mrichter/AliRoot.git] / PWG / DevNanoAOD / AliNanoAODReplicator.cxx

/**************************************************************************
* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
*                                                                        *
* Author: The ALICE Off-line Project.                                    *
* Contributors are mentioned in the code where appropriate.              *
*                                                                        *
* Permission to use, copy, modify and distribute this software and its   *
* documentation strictly for non-commercial purposes is hereby granted   *
* without fee, provided that the above copyright notice appears in all   *
* copies and that both the copyright notice and this permission notice   *
* appear in the supporting documentation. The authors make no claims     *
* about the suitability of this software for any purpose. It is          *
* provided "as is" without express or implied warranty.                  *
**************************************************************************/

// $Id: AliNanoAODReplicator.cxx 56492 2012-05-15 18:42:47Z pcrochet $

// Implementation of a branch replicator 
// to produce nanoAODs.
//
//
// This replicator is in charge of replicating the tracks,vertices,headers
// branches of a standard AOD or ESD file into a nanoAODs 
// (AliAOD.Special.root)
// 
// The class was inspired by AliAODMuonReplicator
// 
// Author: Michele Floris, michele.floris@cern.ch


class AliESDv0;
class AliESDVertex;
class AliAODVertex;
class AliAODRecoDecay;

#include "AliAODDimuon.h"
#include "AliAODEvent.h"
#include "AliAODMCHeader.h"
#include "AliAODMCParticle.h"
#include "AliAODTZERO.h"
#include "AliAODTrack.h"
#include "AliAODVZERO.h"
#include "AliAnalysisCuts.h"
#include "TF1.h"
#include "AliExternalTrackParam.h"
#include "AliESDv0.h"
#include "AliAODv0.h"
#include "AliPIDResponse.h"
#include <iostream>
#include <cassert>
#include "AliESDtrack.h"
#include "TObjArray.h"
#include "AliAnalysisFilter.h"
#include "AliNanoAODTrack.h"

#include <TFile.h>
#include <TDatabasePDG.h>
#include <TString.h>
#include <TList.h>
#include "AliLog.h"
#include "AliVEvent.h"
#include "AliVVertex.h"
#include "AliVTrack.h"
#include "AliVertexerTracks.h"
#include "AliKFVertex.h"
#include "AliESDEvent.h"
#include "AliESDVertex.h"
#include "AliESDtrackCuts.h"
#include "AliAODEvent.h"
#include "AliAnalysisFilter.h"

#include "AliNanoAODReplicator.h"
#include "TH1.h"
#include "TCanvas.h"
#include "AliNanoAODHeader.h"
#include "AliNanoAODCustomSetter.h"

using std::cout;
using std::endl;

ClassImp(AliNanoAODReplicator)

//_____________________________________________________________________________
AliNanoAODReplicator::AliNanoAODReplicator(const char* name, const char* title,
					   const char * varlist,
					   AliAnalysisCuts* trackCut,
					   Int_t mcMode
					     ) :
  AliAODBranchReplicator(name,title), 

  fTrackCut(trackCut), fTracks(0x0), fHeader(0x0), fNTracksVariables(0), // FIXME: Start using cuts, and check if fNTracksVariables is needed
  fVertices(0x0), 
  fList(0x0),
  fMCParticles(0x0),
  fMCHeader(0x0),
  fMCMode(mcMode),
  fLabelMap(),
  fParticleSelected(),
  fVarList(varlist),
  fVarListHeader(""),// FIXME: this should be set to a meaningful value: add an arg to the constructor
  fCustomSetter(0)
{
  // default ctor
  AliNanoAODTrackMapping * tm =new AliNanoAODTrackMapping(fVarList);
  fNTracksVariables = tm->GetSize();  
  //  tm->Print();
}

//_____________________________________________________________________________
AliNanoAODReplicator::~AliNanoAODReplicator()
{
  // dtor
  delete fTrackCut;
  delete fList;
}

//_____________________________________________________________________________
void AliNanoAODReplicator::SelectParticle(Int_t i)
{
  // taking the absolute values here, need to take care 
  // of negative daughter and mother
  // IDs when setting!
  
  if (!IsParticleSelected(TMath::Abs(i)))
  {
    fParticleSelected.Add(TMath::Abs(i),1);    
  }
}

//_____________________________________________________________________________
Bool_t AliNanoAODReplicator::IsParticleSelected(Int_t i)  
{
  // taking the absolute values here, need to take 
  // care with negative daughter and mother
  // IDs when setting!
  return (fParticleSelected.GetValue(TMath::Abs(i))==1);
}


//_____________________________________________________________________________
void AliNanoAODReplicator::CreateLabelMap(const AliAODEvent& source)
{  
  //
  // this should be called once all selections are done 
  // This method associates to the original index of the mc particle
  // (i) the new one (j). J runs only over particles which are
  // actually kept.
  //
  
  fLabelMap.Delete();
  
  TClonesArray* mcParticles = static_cast<TClonesArray*>(source.FindListObject(AliAODMCParticle::StdBranchName()));
  
  Int_t j(0);
  Int_t i(0); // We need i, we cannot rely on part->GetLabel, because some of the original mc particles are not kept in the stack, apparently
  
  TIter next(mcParticles);

  while ( next() )
  {
    if (IsParticleSelected(i))
    {
      fLabelMap.Add(i,j++);
      //      std::cout << i <<  "->" << j-1 << std::endl;
    }
    ++i;
  }  


}

//_____________________________________________________________________________
Int_t AliNanoAODReplicator::GetNewLabel(Int_t i) 
{
  // Gets the label from the new created Map
  // Call CreatLabelMap before
  // otherwise only 0 returned
  return fLabelMap.GetValue(TMath::Abs(i));
}

//_____________________________________________________________________________
void AliNanoAODReplicator::FilterMC(const AliAODEvent& source)
{
  // Filter MC information


  AliAODMCHeader* mcHeader(0x0);
  TClonesArray* mcParticles(0x0);
  
  fParticleSelected.Delete();

  //  std::cout << "MC Mode: " << fMCMode << ", Tracks " << fTracks->GetEntries() << std::endl;
  
  if ( fMCMode>=2 && !fTracks->GetEntries() ) {
    return;
  }
  // for fMCMode==1 we only copy MC information for events where there's at least one muon track
    
  mcHeader = static_cast<AliAODMCHeader*>(source.FindListObject(AliAODMCHeader::StdBranchName()));
  
  if ( mcHeader ) 
    {
      *fMCHeader = *mcHeader;
    }
  
  
  mcParticles = static_cast<TClonesArray*>(source.FindListObject(AliAODMCParticle::StdBranchName()));
  
  if ( mcParticles && fMCMode>=2 )
    {
      // keep all primaries
      TIter nextPart(mcParticles);
      static Int_t iev = -1; // FIXME: remove this (debug)
      iev ++;
      AliAODMCParticle * prim = 0;
      Int_t iprim = 0;  // We need iprim, we cannot rely on part->GetLabel, because some of the original mc particles are not kept in the stack, apparently
      // also select all charged primaries 
      while ((prim = (AliAODMCParticle*) nextPart())) {
	if(prim->IsPhysicalPrimary() && prim->Charge()) SelectParticle(iprim);
	// FIXME DEBUG
	if(iev == 2009) {
	  // std::cout << "IEV " << iev << std::endl;
	  // std::cout << " PART " << iprim << " " << prim->IsPhysicalPrimary() <<","<<prim->Charge() << "=" << IsParticleSelected(iprim) <<  std::endl;
	  if(iprim == 15) {
	    prim->Print();
	  }
	  
	}
	iprim++;
      } 

      // loop on (kept) tracks to find their ancestors
      TIter nextTRACK(fTracks);
      AliNanoAODTrack* track;
    
      while ( ( track = static_cast<AliNanoAODTrack*>(nextTRACK()) ) )
	{
	  Int_t label = TMath::Abs(track->GetLabel()); 
      
	  while ( label >= 0 ) 
	    {
	      SelectParticle(label);
	      AliAODMCParticle* mother = static_cast<AliAODMCParticle*>(mcParticles->UncheckedAt(label));
	      if (!mother)
		{
		  AliError(Form("Got a null mother ! Check that ! (label %d",label)); // FIXME: I think this error is not needed
		  label = -1;
		}
	      else
		{
		  label = mother->GetMother();// do not only keep particles which created a track, but all their mothers
		}
	    }
	}
    
      CreateLabelMap(source);
    
      // Actual filtering and label remapping (shamelessly taken for the implementation of AliAODHandler::StoreMCParticles)
      TIter nextMC(mcParticles);
      AliAODMCParticle* p;
      Int_t nmc(0);  // We need nmc, we cannot rely on part->GetLabel, because some of the original mc particles are not kept in the stack, apparently
      Int_t nmcout(0);
    
      while ( ( p = static_cast<AliAODMCParticle*>(nextMC()) ) )
	{
	  AliAODMCParticle c(*p);

	  if ( IsParticleSelected(nmc) )
	    {
	      // 
	      Int_t d0 =  p->GetDaughter(0);
	      Int_t d1 =  p->GetDaughter(1);
	      Int_t m =   p->GetMother();
        
	      // other than for the track labels, negative values mean
	      // no daughter/mother so preserve it
        
	      if(d0<0 && d1<0)
		{
		  // no first daughter -> no second daughter
		  // nothing to be done
		  // second condition not needed just for sanity check at the end
		  c.SetDaughter(0,d0);
		  c.SetDaughter(1,d1);
		} else if(d1 < 0 && d0 >= 0) 
		{
		  // Only one daughter
		  // second condition not needed just for sanity check at the end
		  if(IsParticleSelected(d0))
		    {
		      c.SetDaughter(0,GetNewLabel(d0));
		    } else 
		    {
		      c.SetDaughter(0,-1);
		    }
		  c.SetDaughter(1,d1);
		}
	      else if (d0 > 0 && d1 > 0 )
		{
		  // we have two or more daughters loop on the stack to see if they are
		  // selected
		  Int_t d0tmp = -1;
		  Int_t d1tmp = -1;
		  for (int id = d0; id<=d1;++id)
		    {
		      if (IsParticleSelected(id))
			{
			  if(d0tmp==-1)
			    {
			      // first time
			      d0tmp = GetNewLabel(id);
			      d1tmp = d0tmp; // this is to have the same schema as on the stack i.e. with one daugther d0 and d1 are the same 
			    }
			  else d1tmp = GetNewLabel(id);
			}
		    }
		  c.SetDaughter(0,d0tmp);
		  c.SetDaughter(1,d1tmp);
		} else 
		{
		  AliFatal(Form("Unxpected indices %d %d",d0,d1));
		}
        
	      if ( m < 0 )
		{
		  c.SetMother(m);
		} else 
		{
		  if (IsParticleSelected(m)) 
		    {
		      c.SetMother(GetNewLabel(m));              
		    }
		  // else // FIXME: re-enable this checj. Sometimes it gets here. Still to be understood why
		  //   {
		  //     //		      AliError(Form("PROBLEM Mother not selected %d", m));              
		  //   }
		}
        
	      new ((*fMCParticles)[nmcout++]) AliAODMCParticle(c);
	    }
      
	  ++nmc;        
	} //closes loop over MC particles
    
      // now remap the tracks...
    
      TIter nextTrack(fTracks);
      AliNanoAODTrack* t;
      //      std::cout << "Remapping tracks" << std::endl;
    
      while ( ( t = dynamic_cast<AliNanoAODTrack*>(nextTrack()) ) )
	{
	  
	  t->SetLabel(GetNewLabel(t->GetLabel()));
	}
    
    } // closes fMCMode == 1
  else if ( mcParticles ) 
    {
      // simple copy of input MC particles to ouput MC particles
      TIter nextMC(mcParticles);
      AliAODMCParticle* p;
      Int_t nmcout(0);
    
      while ( ( p = static_cast<AliAODMCParticle*>(nextMC()) ) )
	{
	  new ((*fMCParticles)[nmcout++]) AliAODMCParticle(*p);
	}
    }
  
  AliDebug(1,Form("input mc %d output mc %d",
                  mcParticles ? mcParticles->GetEntries() : 0,
                  fMCParticles ? fMCParticles->GetEntries() : 0));
  
  Printf("input mc %d output mc %d",
                  mcParticles ? mcParticles->GetEntries() : 0,
                  fMCParticles ? fMCParticles->GetEntries() : 0);
  

}

// //_____________________________________________________________________________
TList* AliNanoAODReplicator::GetList() const
{
  // return (and build if not already done) our internal list of managed objects
  
  if (!fList)
    {
      fList = new TList;
      fList->SetOwner(kTRUE);

      fTracks = new TClonesArray("AliNanoAODTrack");      
      fTracks->SetName("tracks"); // TODO: consider the possibility to use a different name to distinguish in AliAODEvent
      fList->Add(fTracks);    

      fHeader = new AliNanoAODHeader(3);// TODO: to be customized
      fHeader->SetName("header"); // TODO: consider the possibility to use a different name to distinguish in AliAODEvent
      fList->Add(fHeader);    


      fVertices = new TClonesArray("AliAODVertex",2);
      fVertices->SetName("vertices");    
    
        
      fList->Add(fVertices);
    
      if ( fMCMode > 0 )
	{
	  fMCHeader = new AliAODMCHeader;    
	  fMCParticles = new TClonesArray("AliAODMCParticle",1000);
	  fMCParticles->SetName(AliAODMCParticle::StdBranchName());
	  fList->Add(fMCHeader);
	  fList->Add(fMCParticles);
	}
    }
  return fList;
}

//_____________________________________________________________________________
void AliNanoAODReplicator::ReplicateAndFilter(const AliAODEvent& source)
//void AliNanoAODReplicator::ReplicateAndFilter(AliAODEvent *source)
{
  // Replicate (and filter if filters are there) the relevant parts we're interested in AODEvent
  

  // assert(fTracks!=0x0);
  
  //*fTZERO = *(source.GetTZEROData());
  
  
  fTracks->Clear("C");			
  assert(fVertices!=0x0);
  fVertices->Clear("C");
  fMCHeader->Reset();
  fMCParticles->Clear("C");
  
  Int_t ntracks(0);
  Int_t input(0);

  AliAODVertex *vtx = source.GetPrimaryVertex();

  // TODO: implement header!
  //  *fHeader = *source.GetHeader();
  if(fCustomSetter){
    // Set custom variables in the header if the callback is set
    fCustomSetter->SetNanoAODHeader(&source, fHeader);
  }


  //  Int_t nindices=0;
  const Int_t entries = source.GetNumberOfTracks();

  Double_t pos[3],cov[6];
  vtx->GetXYZ(pos);
  vtx->GetCovarianceMatrix(cov);
  
  if(entries<=0) return;

  if(vtx->GetNContributors()<1) {
    // SPD vertex cut
    vtx =source.GetPrimaryVertexSPD();    
    if(vtx->GetNContributors()<1) { // FIXME: Keep this cut?
      Info("AliNanoAODReplicator","No good vertex, skip event");
      return; // NO GOOD VERTEX, SKIP EVENT 
    }
  }
  
  
  for(Int_t j=0; j<entries; j++){
    
    AliVTrack *track = (AliVTrack*)source.GetTrack(j);
    
    AliAODTrack *aodtrack =(AliAODTrack*)track;// FIXME DYNAMIC CAST?
    if(!fTrackCut->IsSelected(aodtrack)) continue;

    AliNanoAODTrack * special = new  AliNanoAODTrack (aodtrack, fVarList);
    
    if(fCustomSetter) fCustomSetter->SetNanoAODTrack(aodtrack, special);
    (*fTracks)[ntracks++] = special;
    //new((*fTracks)[ntrac\ks++])
  }  
  //----------------------------------------------------------
  
  TIter nextV(source.GetVertices());
  AliAODVertex* v;
  Int_t nvertices(0);
  
  while ( ( v = static_cast<AliAODVertex*>(nextV()) ) )
    {
      AliAODVertex* tmp = v->CloneWithoutRefs();
      AliAODVertex* copiedVertex = new((*fVertices)[nvertices++]) AliAODVertex(*tmp);
      
      // to insure the main vertex retains the ncontributors information
      // (which is otherwise computed dynamically from
      // references to tracks, which we do not keep in muon aods...)
      // we set it here
      
      copiedVertex->SetNContributors(v->GetNContributors()); 
      
      //  fVertices->Delete();
      //	  delete copiedVertex;
      delete tmp;
    }
  
  
  AliDebug(1,Form("input mu tracks=%d tracks=%d vertices=%d",
                  input,fTracks->GetEntries(),fVertices->GetEntries())); 
  
  
  // Finally, deal with MC information, if needed
  
  if ( fMCMode > 0 ) {
    FilterMC(source);      
  }
  

}


//-----------------------------------------------------------------------------

//----------------------------------------------------------------------------
	

//-----------------------------------------------------------------------------

// AliAODVertex* AliNanoAODReplicator::PrimaryVertex(const TObjArray *trkArray,
// 						   AliAODEvent &event) const
// {
//   // Returns primary vertex to be used for this candidate
//   //AliCodeTimerAuto("",0);

//   AliESDVertex *vertexESD = 0;
//   AliAODVertex *vertexAOD = 0;


//   if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) { 
//     // primary vertex from the input event
    
//     vertexESD = new AliESDVertex(*fV1);

//   } else {
//     // primary vertex specific to this candidate

//     Int_t nTrks = trkArray->GetEntriesFast();
//     AliVertexerTracks *vertexer = new AliVertexerTracks(event.GetMagneticField());

//     if(fRecoPrimVtxSkippingTrks) { 
//       // recalculating the vertex
      
//       if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
// 	Float_t diamondcovxy[3];
// 	event.GetDiamondCovXY(diamondcovxy);
// 	Double_t pos[3]={event.GetDiamondX(),event.GetDiamondY(),0.};
// 	Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
// 	AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
// 	vertexer->SetVtxStart(diamond);
// 	delete diamond; diamond=NULL;
// 	if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter")) 
// 	  vertexer->SetOnlyFitter();
//       }
//       Int_t skipped[1000];
//       Int_t nTrksToSkip=0,id;
//       AliExternalTrackParam *t = 0;
//       for(Int_t i=0; i<nTrks; i++) {
// 	t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
// 	id = (Int_t)t->GetID();
// 	if(id<0) continue;
// 	skipped[nTrksToSkip++] = id;
//       }
//       // TEMPORARY FIX
//       // For AOD, skip also tracks without covariance matrix
//       if(fInputAOD) {
// 	Double_t covtest[21];
// 	for(Int_t j=0; j<event.GetNumberOfTracks(); j++) {
// 	  AliVTrack *vtrack = (AliVTrack*)event.GetTrack(j);
// 	  if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
// 	    id = (Int_t)vtrack->GetID();
// 	    if(id<0) continue;
// 	    skipped[nTrksToSkip++] = id;
// 	  }
// 	}
//       }
//       for(Int_t ijk=nTrksToSkip; ijk<1000; ijk++) skipped[ijk]=-1;
//       //
//       vertexer->SetSkipTracks(nTrksToSkip,skipped);
//       vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event); 
      
//     } else if(fRmTrksFromPrimVtx && nTrks>0) { 
//       // removing the prongs tracks
      
//       TObjArray rmArray(nTrks);
//       UShort_t *rmId = new UShort_t[nTrks];
//       AliESDtrack *esdTrack = 0;
//       AliESDtrack *t = 0;
//       for(Int_t i=0; i<nTrks; i++) {
// 	t = (AliESDtrack*)trkArray->UncheckedAt(i);
// 	esdTrack = new AliESDtrack(*t);
// 	rmArray.AddLast(esdTrack);
// 	if(esdTrack->GetID()>=0) {
// 	  rmId[i]=(UShort_t)esdTrack->GetID();
// 	} else {
// 	  rmId[i]=9999;
// 	}
//       }
//       Float_t diamondxy[2]={event.GetDiamondX(),event.GetDiamondY()};
//       vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
//       delete [] rmId; rmId=NULL;
//       rmArray.Delete();
      
//     }

//     if(!vertexESD) return vertexAOD;
//     if(vertexESD->GetNContributors()<=0) { 
//       //AliDebug(2,"vertexing failed"); 
//       delete vertexESD; vertexESD=NULL;
//       return vertexAOD;
//     }

//     delete vertexer; vertexer=NULL;

//   }

//   // convert to AliAODVertex
//   Double_t pos[3],cov[6],chi2perNDF;
//   vertexESD->GetXYZ(pos); // position
//   vertexESD->GetCovMatrix(cov); //covariance matrix
//   chi2perNDF = vertexESD->GetChi2toNDF();
//   delete vertexESD; vertexESD=NULL;

//   vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);

//   return vertexAOD;
// }

//_____________________________________________________________________________


// //---------------------------------------------------------------------------

void AliNanoAODReplicator::Terminate(){

}
Commit	Line	Data
778636d4	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
	16	// $Id: AliNanoAODReplicator.cxx 56492 2012-05-15 18:42:47Z pcrochet $
	17
	18	// Implementation of a branch replicator
	19	// to produce nanoAODs.
	20	//
	21	//
	22	// This replicator is in charge of replicating the tracks,vertices,headers
	23	// branches of a standard AOD or ESD file into a nanoAODs
	24	// (AliAOD.Special.root)
	25	//
	26	// The class was inspired by AliAODMuonReplicator
	27	//
	28	// Author: Michele Floris, michele.floris@cern.ch
	29
	30
	31	class AliESDv0;
	32	class AliESDVertex;
	33	class AliAODVertex;
	34	class AliAODRecoDecay;
	35
	36	#include "AliAODDimuon.h"
	37	#include "AliAODEvent.h"
	38	#include "AliAODMCHeader.h"
	39	#include "AliAODMCParticle.h"
	40	#include "AliAODTZERO.h"
	41	#include "AliAODTrack.h"
	42	#include "AliAODVZERO.h"
	43	#include "AliAnalysisCuts.h"
	44	#include "TF1.h"
	45	#include "AliExternalTrackParam.h"
	46	#include "AliESDv0.h"
	47	#include "AliAODv0.h"
	48	#include "AliPIDResponse.h"
	49	#include <iostream>
	50	#include <cassert>
	51	#include "AliESDtrack.h"
	52	#include "TObjArray.h"
	53	#include "AliAnalysisFilter.h"
	54	#include "AliNanoAODTrack.h"
	55
	56	#include <TFile.h>
	57	#include <TDatabasePDG.h>
	58	#include <TString.h>
	59	#include <TList.h>
	60	#include "AliLog.h"
	61	#include "AliVEvent.h"
	62	#include "AliVVertex.h"
	63	#include "AliVTrack.h"
	64	#include "AliVertexerTracks.h"
65	#include "AliKFVertex.h"
66	#include "AliESDEvent.h"
67	#include "AliESDVertex.h"
68	#include "AliESDtrackCuts.h"
69	#include "AliAODEvent.h"
70	#include "AliAnalysisFilter.h"
71
72	#include "AliNanoAODReplicator.h"
73	#include "TH1.h"
74	#include "TCanvas.h"
75	#include "AliNanoAODHeader.h"
76	#include "AliNanoAODCustomSetter.h"
77
78	using std::cout;
79	using std::endl;
80
81	ClassImp(AliNanoAODReplicator)
82
83	//_____________________________________________________________________________
84	AliNanoAODReplicator::AliNanoAODReplicator(const char* name, const char* title,
9c35eb62	85	const char * varlist,
	86	AliAnalysisCuts* trackCut,
	87	Int_t mcMode
778636d4	88	) :
	89	AliAODBranchReplicator(name,title),
	90
	91	fTrackCut(trackCut), fTracks(0x0), fHeader(0x0), fNTracksVariables(0), // FIXME: Start using cuts, and check if fNTracksVariables is needed
	92	fVertices(0x0),
	93	fList(0x0),
	94	fMCParticles(0x0),
	95	fMCHeader(0x0),
	96	fMCMode(mcMode),
	97	fLabelMap(),
	98	fParticleSelected(),
	99	fVarList(varlist),
9c35eb62	100	fVarListHeader(""),// FIXME: this should be set to a meaningful value: add an arg to the constructor
778636d4	101	fCustomSetter(0)
	102	{
	103	// default ctor
	104	AliNanoAODTrackMapping * tm =new AliNanoAODTrackMapping(fVarList);
	105	fNTracksVariables = tm->GetSize();
	106	// tm->Print();
	107	}
	108
	109	//_____________________________________________________________________________
	110	AliNanoAODReplicator::~AliNanoAODReplicator()
	111	{
	112	// dtor
	113	delete fTrackCut;
	114	delete fList;
	115	}
	116
	117	//_____________________________________________________________________________
	118	void AliNanoAODReplicator::SelectParticle(Int_t i)
	119	{
	120	// taking the absolute values here, need to take care
	121	// of negative daughter and mother
	122	// IDs when setting!
	123
	124	if (!IsParticleSelected(TMath::Abs(i)))
	125	{
	126	fParticleSelected.Add(TMath::Abs(i),1);
	127	}
	128	}
	129
	130	//_____________________________________________________________________________
	131	Bool_t AliNanoAODReplicator::IsParticleSelected(Int_t i)
	132	{
	133	// taking the absolute values here, need to take
	134	// care with negative daughter and mother
	135	// IDs when setting!
	136	return (fParticleSelected.GetValue(TMath::Abs(i))==1);
	137	}
	138
	139
	140	//_____________________________________________________________________________
	141	void AliNanoAODReplicator::CreateLabelMap(const AliAODEvent& source)
	142	{
	143	//
	144	// this should be called once all selections are done
	145	// This method associates to the original index of the mc particle
	146	// (i) the new one (j). J runs only over particles which are
	147	// actually kept.
	148	//
	149
	150	fLabelMap.Delete();
	151
	152	TClonesArray* mcParticles = static_cast<TClonesArray*>(source.FindListObject(AliAODMCParticle::StdBranchName()));
	153
	154	Int_t j(0);
	155	Int_t i(0); // We need i, we cannot rely on part->GetLabel, because some of the original mc particles are not kept in the stack, apparently
	156
	157	TIter next(mcParticles);
	158
	159	while ( next() )
	160	{
	161	if (IsParticleSelected(i))
	162	{
	163	fLabelMap.Add(i,j++);
	164	// std::cout << i << "->" << j-1 << std::endl;
165	}
166	++i;
167	}
168
169
170	}
171
172	//_____________________________________________________________________________
173	Int_t AliNanoAODReplicator::GetNewLabel(Int_t i)
174	{
175	// Gets the label from the new created Map
176	// Call CreatLabelMap before
177	// otherwise only 0 returned
178	return fLabelMap.GetValue(TMath::Abs(i));
179	}
180
181	//_____________________________________________________________________________
182	void AliNanoAODReplicator::FilterMC(const AliAODEvent& source)
183	{
184	// Filter MC information
185
186
187	AliAODMCHeader* mcHeader(0x0);
188	TClonesArray* mcParticles(0x0);
189
190	fParticleSelected.Delete();
191
192	// std::cout << "MC Mode: " << fMCMode << ", Tracks " << fTracks->GetEntries() << std::endl;
193
194	if ( fMCMode>=2 && !fTracks->GetEntries() ) {
195	return;
196	}
197	// for fMCMode==1 we only copy MC information for events where there's at least one muon track
198
199	mcHeader = static_cast<AliAODMCHeader*>(source.FindListObject(AliAODMCHeader::StdBranchName()));
200
201	if ( mcHeader )
202	{
203	fMCHeader = mcHeader;
204	}
205
206
207	mcParticles = static_cast<TClonesArray*>(source.FindListObject(AliAODMCParticle::StdBranchName()));
208
209	if ( mcParticles && fMCMode>=2 )
210	{
211	// keep all primaries
212	TIter nextPart(mcParticles);
213	static Int_t iev = -1; // FIXME: remove this (debug)
214	iev ++;
215	AliAODMCParticle * prim = 0;
216	Int_t iprim = 0; // We need iprim, we cannot rely on part->GetLabel, because some of the original mc particles are not kept in the stack, apparently
217	// also select all charged primaries
218	while ((prim = (AliAODMCParticle*) nextPart())) {
219	if(prim->IsPhysicalPrimary() && prim->Charge()) SelectParticle(iprim);
220	// FIXME DEBUG
221	if(iev == 2009) {
222	// std::cout << "IEV " << iev << std::endl;
223	// std::cout << " PART " << iprim << " " << prim->IsPhysicalPrimary() <<","<<prim->Charge() << "=" << IsParticleSelected(iprim) << std::endl;
224	if(iprim == 15) {
225	prim->Print();
226	}
227
228	}
229	iprim++;
230	}
231
232	// loop on (kept) tracks to find their ancestors
233	TIter nextTRACK(fTracks);
234	AliNanoAODTrack* track;
235
236	while ( ( track = static_cast<AliNanoAODTrack*>(nextTRACK()) ) )
237	{
238	Int_t label = TMath::Abs(track->GetLabel());
239
240	while ( label >= 0 )
241	{
242	SelectParticle(label);
243	AliAODMCParticle* mother = static_cast<AliAODMCParticle*>(mcParticles->UncheckedAt(label));
244	if (!mother)
245	{
246	AliError(Form("Got a null mother ! Check that ! (label %d",label)); // FIXME: I think this error is not needed
247	label = -1;
248	}
249	else
250	{
251	label = mother->GetMother();// do not only keep particles which created a track, but all their mothers
252	}
253	}
254	}
255
256	CreateLabelMap(source);
257
258	// Actual filtering and label remapping (shamelessly taken for the implementation of AliAODHandler::StoreMCParticles)
259	TIter nextMC(mcParticles);
260	AliAODMCParticle* p;
261	Int_t nmc(0); // We need nmc, we cannot rely on part->GetLabel, because some of the original mc particles are not kept in the stack, apparently
262	Int_t nmcout(0);
263
264	while ( ( p = static_cast<AliAODMCParticle*>(nextMC()) ) )
265	{
266	AliAODMCParticle c(*p);
267
268	if ( IsParticleSelected(nmc) )
269	{
270	//
271	Int_t d0 = p->GetDaughter(0);
272	Int_t d1 = p->GetDaughter(1);
273	Int_t m = p->GetMother();
274
275	// other than for the track labels, negative values mean
276	// no daughter/mother so preserve it
277
278	if(d0<0 && d1<0)
279	{
280	// no first daughter -> no second daughter
281	// nothing to be done
282	// second condition not needed just for sanity check at the end
283	c.SetDaughter(0,d0);
284	c.SetDaughter(1,d1);
285	} else if(d1 < 0 && d0 >= 0)
286	{
287	// Only one daughter
288	// second condition not needed just for sanity check at the end
289	if(IsParticleSelected(d0))
290	{
291	c.SetDaughter(0,GetNewLabel(d0));
292	} else
293	{
294	c.SetDaughter(0,-1);
295	}
296	c.SetDaughter(1,d1);
297	}
298	else if (d0 > 0 && d1 > 0 )
299	{
300	// we have two or more daughters loop on the stack to see if they are
301	// selected
302	Int_t d0tmp = -1;
303	Int_t d1tmp = -1;
304	for (int id = d0; id<=d1;++id)
305	{
306	if (IsParticleSelected(id))
307	{
308	if(d0tmp==-1)
309	{
310	// first time
311	d0tmp = GetNewLabel(id);
312	d1tmp = d0tmp; // this is to have the same schema as on the stack i.e. with one daugther d0 and d1 are the same
313	}
314	else d1tmp = GetNewLabel(id);
315	}
316	}
317	c.SetDaughter(0,d0tmp);
318	c.SetDaughter(1,d1tmp);
319	} else
320	{
321	AliFatal(Form("Unxpected indices %d %d",d0,d1));
322	}
323
324	if ( m < 0 )
325	{
326	c.SetMother(m);
327	} else
328	{
329	if (IsParticleSelected(m))
330	{
331	c.SetMother(GetNewLabel(m));
332	}
333	// else // FIXME: re-enable this checj. Sometimes it gets here. Still to be understood why
334	// {
335	// // AliError(Form("PROBLEM Mother not selected %d", m));
336	// }
337	}
338
339	new ((*fMCParticles)[nmcout++]) AliAODMCParticle(c);
340	}
341
342	++nmc;
343	} //closes loop over MC particles
344
345	// now remap the tracks...
346
347	TIter nextTrack(fTracks);
348	AliNanoAODTrack* t;
349	// std::cout << "Remapping tracks" << std::endl;
350
351	while ( ( t = dynamic_cast<AliNanoAODTrack*>(nextTrack()) ) )
352	{
353
354	t->SetLabel(GetNewLabel(t->GetLabel()));
355	}
356
357	} // closes fMCMode == 1
358	else if ( mcParticles )
359	{
360	// simple copy of input MC particles to ouput MC particles
361	TIter nextMC(mcParticles);
362	AliAODMCParticle* p;
363	Int_t nmcout(0);
364
365	while ( ( p = static_cast<AliAODMCParticle*>(nextMC()) ) )
366	{
367	new ((fMCParticles)[nmcout++]) AliAODMCParticle(p);
368	}
369	}
370
371	AliDebug(1,Form("input mc %d output mc %d",
372	mcParticles ? mcParticles->GetEntries() : 0,
373	fMCParticles ? fMCParticles->GetEntries() : 0));
374
375	Printf("input mc %d output mc %d",
376	mcParticles ? mcParticles->GetEntries() : 0,
377	fMCParticles ? fMCParticles->GetEntries() : 0);
378
379
380	}
381
382	// //_____________________________________________________________________________
383	TList* AliNanoAODReplicator::GetList() const
384	{
385	// return (and build if not already done) our internal list of managed objects
386
387	if (!fList)
388	{
389	fList = new TList;
390	fList->SetOwner(kTRUE);
391
392	fTracks = new TClonesArray("AliNanoAODTrack");
393	fTracks->SetName("tracks"); // TODO: consider the possibility to use a different name to distinguish in AliAODEvent
394	fList->Add(fTracks);
395
396	fHeader = new AliNanoAODHeader(3);// TODO: to be customized
397	fHeader->SetName("header"); // TODO: consider the possibility to use a different name to distinguish in AliAODEvent
398	fList->Add(fHeader);
399
400
401	fVertices = new TClonesArray("AliAODVertex",2);
402	fVertices->SetName("vertices");
403
404
405	fList->Add(fVertices);
406
407	if ( fMCMode > 0 )
408	{
409	fMCHeader = new AliAODMCHeader;
410	fMCParticles = new TClonesArray("AliAODMCParticle",1000);
411	fMCParticles->SetName(AliAODMCParticle::StdBranchName());
412	fList->Add(fMCHeader);
413	fList->Add(fMCParticles);
414	}
415	}
416	return fList;
417	}
418
419	//_____________________________________________________________________________
420	void AliNanoAODReplicator::ReplicateAndFilter(const AliAODEvent& source)
421	//void AliNanoAODReplicator::ReplicateAndFilter(AliAODEvent *source)
422	{
423	// Replicate (and filter if filters are there) the relevant parts we're interested in AODEvent
424
425
426	// assert(fTracks!=0x0);
427
428	//fTZERO = (source.GetTZEROData());
429
430
431
432	fTracks->Clear("C");
433	assert(fVertices!=0x0);
434	fVertices->Clear("C");
435	fMCHeader->Reset();
436	fMCParticles->Clear("C");
437
438	Int_t ntracks(0);
439	Int_t input(0);
440
441	AliAODVertex *vtx = source.GetPrimaryVertex();
442
443	// TODO: implement header!
444	// fHeader = source.GetHeader();
445	if(fCustomSetter){
446	// Set custom variables in the header if the callback is set
447	fCustomSetter->SetNanoAODHeader(&source, fHeader);
448	}
449
450
451	// Int_t nindices=0;
452	const Int_t entries = source.GetNumberOfTracks();
453
454	Double_t pos[3],cov[6];
455	vtx->GetXYZ(pos);
456	vtx->GetCovarianceMatrix(cov);
457
458	if(entries<=0) return;
459
460	if(vtx->GetNContributors()<1) {
461	// SPD vertex cut
462	vtx =source.GetPrimaryVertexSPD();
463	if(vtx->GetNContributors()<1) { // FIXME: Keep this cut?
464	Info("AliNanoAODReplicator","No good vertex, skip event");
465	return; // NO GOOD VERTEX, SKIP EVENT
466	}
467	}
468
778636d4	469
	470	for(Int_t j=0; j<entries; j++){
	471
	472	AliVTrack track = (AliVTrack)source.GetTrack(j);
	473
	474	AliAODTrack aodtrack =(AliAODTrack)track;// FIXME DYNAMIC CAST?
	475	if(!fTrackCut->IsSelected(aodtrack)) continue;
	476
	477	AliNanoAODTrack * special = new AliNanoAODTrack (aodtrack, fVarList);
	478
	479	if(fCustomSetter) fCustomSetter->SetNanoAODTrack(aodtrack, special);
	480	(*fTracks)[ntracks++] = special;
	481	//new((*fTracks)[ntrac\ks++])
	482	}
	483	//----------------------------------------------------------
	484
	485	TIter nextV(source.GetVertices());
	486	AliAODVertex* v;
	487	Int_t nvertices(0);
	488
	489	while ( ( v = static_cast<AliAODVertex*>(nextV()) ) )
	490	{
	491	AliAODVertex* tmp = v->CloneWithoutRefs();
	492	AliAODVertex* copiedVertex = new((fVertices)[nvertices++]) AliAODVertex(tmp);
	493
	494	// to insure the main vertex retains the ncontributors information
	495	// (which is otherwise computed dynamically from
	496	// references to tracks, which we do not keep in muon aods...)
	497	// we set it here
	498
	499	copiedVertex->SetNContributors(v->GetNContributors());
	500
	501	// fVertices->Delete();
	502	// delete copiedVertex;
	503	delete tmp;
	504	}
	505
	506
	507	AliDebug(1,Form("input mu tracks=%d tracks=%d vertices=%d",
	508	input,fTracks->GetEntries(),fVertices->GetEntries()));
	509
	510
	511	// Finally, deal with MC information, if needed
	512
	513	if ( fMCMode > 0 ) {
	514	FilterMC(source);
	515	}
	516
	517
	518	}
	519
	520
	521
	522	//-----------------------------------------------------------------------------
	523
	524	//----------------------------------------------------------------------------
	525
	526
	527	//-----------------------------------------------------------------------------
	528
	529	// AliAODVertex* AliNanoAODReplicator::PrimaryVertex(const TObjArray *trkArray,
	530	// AliAODEvent &event) const
	531	// {
	532	// // Returns primary vertex to be used for this candidate
533	// //AliCodeTimerAuto("",0);
534
535	// AliESDVertex *vertexESD = 0;
536	// AliAODVertex *vertexAOD = 0;
537
538
539	// if(!fRecoPrimVtxSkippingTrks && !fRmTrksFromPrimVtx) {
540	// // primary vertex from the input event
541
542	// vertexESD = new AliESDVertex(*fV1);
543
544	// } else {
545	// // primary vertex specific to this candidate
546
547	// Int_t nTrks = trkArray->GetEntriesFast();
548	// AliVertexerTracks *vertexer = new AliVertexerTracks(event.GetMagneticField());
549
550	// if(fRecoPrimVtxSkippingTrks) {
551	// // recalculating the vertex
552
553	// if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraint")) {
554	// Float_t diamondcovxy[3];
555	// event.GetDiamondCovXY(diamondcovxy);
556	// Double_t pos[3]={event.GetDiamondX(),event.GetDiamondY(),0.};
557	// Double_t cov[6]={diamondcovxy[0],diamondcovxy[1],diamondcovxy[2],0.,0.,10.*10.};
558	// AliESDVertex *diamond = new AliESDVertex(pos,cov,1.,1);
559	// vertexer->SetVtxStart(diamond);
560	// delete diamond; diamond=NULL;
561	// if(strstr(fV1->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
562	// vertexer->SetOnlyFitter();
563	// }
564	// Int_t skipped[1000];
565	// Int_t nTrksToSkip=0,id;
566	// AliExternalTrackParam *t = 0;
567	// for(Int_t i=0; i<nTrks; i++) {
568	// t = (AliExternalTrackParam*)trkArray->UncheckedAt(i);
569	// id = (Int_t)t->GetID();
570	// if(id<0) continue;
571	// skipped[nTrksToSkip++] = id;
572	// }
573	// // TEMPORARY FIX
574	// // For AOD, skip also tracks without covariance matrix
575	// if(fInputAOD) {
576	// Double_t covtest[21];
577	// for(Int_t j=0; j<event.GetNumberOfTracks(); j++) {
578	// AliVTrack vtrack = (AliVTrack)event.GetTrack(j);
579	// if(!vtrack->GetCovarianceXYZPxPyPz(covtest)) {
580	// id = (Int_t)vtrack->GetID();
581	// if(id<0) continue;
582	// skipped[nTrksToSkip++] = id;
583	// }
584	// }
585	// }
586	// for(Int_t ijk=nTrksToSkip; ijk<1000; ijk++) skipped[ijk]=-1;
587	// //
588	// vertexer->SetSkipTracks(nTrksToSkip,skipped);
589	// vertexESD = (AliESDVertex*)vertexer->FindPrimaryVertex(event);
590
591	// } else if(fRmTrksFromPrimVtx && nTrks>0) {
592	// // removing the prongs tracks
593
594	// TObjArray rmArray(nTrks);
595	// UShort_t *rmId = new UShort_t[nTrks];
596	// AliESDtrack *esdTrack = 0;
597	// AliESDtrack *t = 0;
598	// for(Int_t i=0; i<nTrks; i++) {
599	// t = (AliESDtrack*)trkArray->UncheckedAt(i);
600	// esdTrack = new AliESDtrack(*t);
601	// rmArray.AddLast(esdTrack);
602	// if(esdTrack->GetID()>=0) {
603	// rmId[i]=(UShort_t)esdTrack->GetID();
604	// } else {
605	// rmId[i]=9999;
606	// }
607	// }
608	// Float_t diamondxy[2]={event.GetDiamondX(),event.GetDiamondY()};
609	// vertexESD = vertexer->RemoveTracksFromVertex(fV1,&rmArray,rmId,diamondxy);
610	// delete [] rmId; rmId=NULL;
611	// rmArray.Delete();
612
613	// }
614
615	// if(!vertexESD) return vertexAOD;
616	// if(vertexESD->GetNContributors()<=0) {
617	// //AliDebug(2,"vertexing failed");
618	// delete vertexESD; vertexESD=NULL;
619	// return vertexAOD;
620	// }
621
622	// delete vertexer; vertexer=NULL;
623
624	// }
625
626	// // convert to AliAODVertex
627	// Double_t pos[3],cov[6],chi2perNDF;
628	// vertexESD->GetXYZ(pos); // position
629	// vertexESD->GetCovMatrix(cov); //covariance matrix
630	// chi2perNDF = vertexESD->GetChi2toNDF();
631	// delete vertexESD; vertexESD=NULL;
632
633	// vertexAOD = new AliAODVertex(pos,cov,chi2perNDF);
634
635	// return vertexAOD;
636	// }
637
638	//_____________________________________________________________________________
639
640
641
642	// //---------------------------------------------------------------------------
643
644	void AliNanoAODReplicator::Terminate(){
645
646	}