PWG3/muon/AliAODMuonReplicator.cxx

   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$
  17
  18 //
  19 // Implementation of a branch replicator
  20 // to produce slim muon and dimuon aods.
  21 //
  22 // This replicator is in charge of replicating the tracks,vertices and dimuons
  23 // branches of the standard AOD into muon AODs (AliAOD.Muons.root and
  24 // AliAOD.Dimuons.root)
  25 //
  26 // The tracks are filtered so that only muon tracks (and only muon tracks
  27 // that pass the trackCut if present) make it to the output aods
  28 //
  29 // The vertices are filtered so that only the primary vertices make it
  30 // to the output aods.
  31 //
  32 // The dimuons are recreated here, according to the set of tracks
  33 // that pass the trackCut (that set may be the same as the input set,
  34 // but to be 100% safe, we always recreate the dimuons).
  35 //
  36 // Author: L. Aphecetche (Subatech)
  37
  38 #include "AliAODMuonReplicator.h"
  39 #include "AliAnalysisCuts.h"
  40 #include "AliAODEvent.h"
  41 #include "AliAODTrack.h"
  42 #include "AliAODDimuon.h"
  43 #include "AliAODMCHeader.h"
  44 #include "AliAODMCParticle.h"
  45 #include <cassert>
  46
  47 ClassImp(AliAODMuonReplicator)
  48
  49 //_____________________________________________________________________________
  50 AliAODMuonReplicator::AliAODMuonReplicator(const char* name, const char* title,
  51                                            AliAnalysisCuts* trackCut,
  52                                            AliAnalysisCuts* vertexCut,
  53                                            Int_t mcMode)
  54 :AliAODBranchReplicator(name,title),
  55 fTrackCut(trackCut), fTracks(0x0),
  56 fVertexCut(vertexCut), fVertices(0x0),
  57 fDimuons(0x0),
  58 fList(0x0),
  59 fMCParticles(0x0),
  60 fMCHeader(0x0),
  61 fMCMode(mcMode),
  62 fLabelMap(),
  63 fParticleSelected()
  64 {
  65   // default ctor
  66 }
  67
  68 //_____________________________________________________________________________
  69 AliAODMuonReplicator::~AliAODMuonReplicator()
  70 {
  71   // dtor
  72   delete fTrackCut;
  73   delete fVertexCut;
  74   delete fList;
  75 }
  76
  77 //_____________________________________________________________________________
  78 void AliAODMuonReplicator::SelectParticle(Int_t i)
  79 {
  80   // taking the absolute values here, need to take care
  81   // of negative daughter and mother
  82   // IDs when setting!
  83
  84   if (!IsParticleSelected(TMath::Abs(i)))
  85   {
  86     fParticleSelected.Add(TMath::Abs(i),1);
  87   }
  88 }
  89
  90 //_____________________________________________________________________________
  91 Bool_t AliAODMuonReplicator::IsParticleSelected(Int_t i)
  92 {
  93   // taking the absolute values here, need to take
  94   // care with negative daughter and mother
  95   // IDs when setting!
  96   return (fParticleSelected.GetValue(TMath::Abs(i))==1);
  97 }
  98
  99
 100 //_____________________________________________________________________________
 101 void AliAODMuonReplicator::CreateLabelMap(const AliAODEvent& source)
 102 {
 103   //
 104   // this should be called once all selections are done
 105   //
 106
 107   fLabelMap.Delete();
 108
 109   TClonesArray* mcParticles = static_cast<TClonesArray*>(source.FindListObject(AliAODMCParticle::StdBranchName()));
 110
 111   Int_t i(0);
 112   Int_t j(0);
 113
 114   TIter next(mcParticles);
 115
 116   while ( next() )
 117   {
 118     if (IsParticleSelected(i))
 119     {
 120       fLabelMap.Add(i,j++);
 121     }
 122     ++i;
 123   }
 124 }
 125
 126 //_____________________________________________________________________________
 127 Int_t AliAODMuonReplicator::GetNewLabel(Int_t i)
 128 {
 129   // Gets the label from the new created Map
 130   // Call CreatLabelMap before
 131   // otherwise only 0 returned
 132   return fLabelMap.GetValue(TMath::Abs(i));
 133 }
 134
 135 //_____________________________________________________________________________
 136 void AliAODMuonReplicator::FilterMC(const AliAODEvent& source)
 137 {
 138   // Filter MC information
 139
 140   fMCHeader->Reset();
 141   fMCParticles->Clear("C");
 142
 143   AliAODMCHeader* mcHeader(0x0);
 144   TClonesArray* mcParticles(0x0);
 145
 146   fParticleSelected.Delete();
 147
 148   if ( !fTracks->GetEntries() ) return;
 149   // only copy MC information for events where there's at least one muon track
 150
 151   mcHeader = static_cast<AliAODMCHeader*>(source.FindListObject(AliAODMCHeader::StdBranchName()));
 152
 153   if ( mcHeader )
 154   {
 155     *fMCHeader = *mcHeader;
 156   }
 157
 158   mcParticles = static_cast<TClonesArray*>(source.FindListObject(AliAODMCParticle::StdBranchName()));
 159
 160   if ( mcParticles && fMCMode>=2 )
 161   {
 162     // loop on (kept) muon tracks to find their ancestors
 163     TIter nextMT(fTracks);
 164     AliAODTrack* mt;
 165
 166     while ( ( mt = static_cast<AliAODTrack*>(nextMT()) ) )
 167     {
 168       Int_t label = mt->GetLabel();
 169
 170       while ( label >= 0 )
 171       {
 172         SelectParticle(label);
 173         AliAODMCParticle* mother = static_cast<AliAODMCParticle*>(mcParticles->UncheckedAt(label));
 174         if (!mother)
 175         {
 176           AliError("Got a null mother ! Check that !");
 177           label = -1;
 178         }
 179         else
 180         {
 181           label = mother->GetMother();
 182         }
 183       }
 184     }
 185
 186     CreateLabelMap(source);
 187
 188     // Actual filtering and label remapping (shamelessly taken for the implementation of AliAODHandler::StoreMCParticles)
 189     TIter nextMC(mcParticles);
 190     AliAODMCParticle* p;
 191     Int_t nmc(0);
 192     Int_t nmcout(0);
 193
 194     while ( ( p = static_cast<AliAODMCParticle*>(nextMC()) ) )
 195     {
 196       AliAODMCParticle c(*p);
 197
 198       if ( IsParticleSelected(nmc) )
 199       {
 200         //
 201         Int_t d0 =  p->GetDaughter(0);
 202         Int_t d1 =  p->GetDaughter(1);
 203         Int_t m =   p->GetMother();
 204
 205         // other than for the track labels, negative values mean
 206         // no daughter/mother so preserve it
 207
 208         if(d0<0 && d1<0)
 209         {
 210           // no first daughter -> no second daughter
 211           // nothing to be done
 212           // second condition not needed just for sanity check at the end
 213           c.SetDaughter(0,d0);
 214           c.SetDaughter(1,d1);
 215         } else if(d1 < 0 && d0 >= 0)
 216         {
 217           // Only one daughter
 218           // second condition not needed just for sanity check at the end
 219           if(IsParticleSelected(d0))
 220           {
 221             c.SetDaughter(0,GetNewLabel(d0));
 222           } else
 223           {
 224             c.SetDaughter(0,-1);
 225           }
 226           c.SetDaughter(1,d1);
 227         }
 228         else if (d0 > 0 && d1 > 0 )
 229         {
 230           // we have two or more daughters loop on the stack to see if they are
 231           // selected
 232           Int_t d0tmp = -1;
 233           Int_t d1tmp = -1;
 234           for (int id = d0; id<=d1;++id)
 235           {
 236             if (IsParticleSelected(id))
 237             {
 238               if(d0tmp==-1)
 239               {
 240                 // first time
 241                 d0tmp = GetNewLabel(id);
 242                 d1tmp = d0tmp; // this is to have the same schema as on the stack i.e. with one daugther d0 and d1 are the same
 243               }
 244               else d1tmp = GetNewLabel(id);
 245             }
 246           }
 247           c.SetDaughter(0,d0tmp);
 248           c.SetDaughter(1,d1tmp);
 249         } else
 250         {
 251           AliError(Form("Unxpected indices %d %d",d0,d1));
 252         }
 253
 254         if ( m < 0 )
 255         {
 256           c.SetMother(m);
 257         } else
 258         {
 259           if (IsParticleSelected(m))
 260           {
 261             c.SetMother(GetNewLabel(m));
 262           }
 263           else
 264           {
 265             AliError(Form("PROBLEM Mother not selected %d", m));
 266           }
 267         }
 268
 269         new ((*fMCParticles)[nmcout++]) AliAODMCParticle(c);
 270       }
 271
 272       ++nmc;
 273     }
 274
 275     // now remap the tracks...
 276
 277     TIter nextTrack(fTracks);
 278     AliAODTrack* t;
 279
 280     while ( ( t = static_cast<AliAODTrack*>(nextTrack()) ) )
 281     {
 282       t->SetLabel(GetNewLabel(t->GetLabel()));
 283     }
 284
 285   }
 286   else if ( mcParticles )
 287   {
 288     // simple copy of input MC particles to ouput MC particles
 289     TIter nextMC(mcParticles);
 290     AliAODMCParticle* p;
 291     Int_t nmcout(0);
 292
 293     while ( ( p = static_cast<AliAODMCParticle*>(nextMC()) ) )
 294     {
 295       new ((*fMCParticles)[nmcout++]) AliAODMCParticle(*p);
 296     }
 297   }
 298
 299   AliDebug(1,Form("input mc %d output mc %d",
 300                   mcParticles ? mcParticles->GetEntries() : 0,
 301                   fMCParticles ? fMCParticles->GetEntries() : 0));
 302
 303 }
 304
 305 //_____________________________________________________________________________
 306 TList* AliAODMuonReplicator::GetList() const
 307 {
 308   // return (and build if not already done) our internal list of managed objects
 309   if (!fList)
 310   {
 311     fTracks = new TClonesArray("AliAODTrack",30);
 312                 fTracks->SetName("tracks");
 313
 314     fVertices = new TClonesArray("AliAODVertex",2);
 315                 fVertices->SetName("vertices");
 316
 317     fDimuons = new TClonesArray("AliAODDimuon",2);
 318     fDimuons->SetName("dimuons");
 319
 320     fList = new TList;
 321     fList->SetOwner(kTRUE);
 322
 323     fList->Add(fTracks);
 324     fList->Add(fVertices);
 325     fList->Add(fDimuons);
 326
 327     if ( fMCMode > 0 )
 328     {
 329       fMCHeader = new AliAODMCHeader;
 330       fMCParticles = new TClonesArray("AliAODMCParticle",1000);
 331       fMCParticles->SetName(AliAODMCParticle::StdBranchName());
 332       fList->Add(fMCHeader);
 333       fList->Add(fMCParticles);
 334     }
 335   }
 336   return fList;
 337 }
 338
 339 //_____________________________________________________________________________
 340 void AliAODMuonReplicator::ReplicateAndFilter(const AliAODEvent& source)
 341 {
 342   // Replicate (and filter if filters are there) the relevant parts we're interested in AODEvent
 343
 344   assert(fTracks!=0x0);
 345   fTracks->Clear("C");
 346   TIter next(source.GetTracks());
 347   AliAODTrack* t;
 348   Int_t ntracks(0);
 349   Int_t input(0);
 350
 351   while ( ( t = static_cast<AliAODTrack*>(next()) ) )
 352   {
 353     if ( t->IsMuonTrack() )
 354     {
 355       ++input;
 356     }
 357
 358     if ( !fTrackCut || fTrackCut->IsSelected(t) )
 359     {
 360       new((*fTracks)[ntracks++]) AliAODTrack(*t);
 361     }
 362   }
 363
 364   assert(fVertices!=0x0);
 365   fVertices->Clear("C");
 366   TIter nextV(source.GetVertices());
 367   AliAODVertex* v;
 368   Int_t nvertices(0);
 369
 370   while ( ( v = static_cast<AliAODVertex*>(nextV()) ) )
 371   {
 372     if ( !fVertexCut || fVertexCut->IsSelected(v) )
 373     {
 374       AliAODVertex* tmp = v->CloneWithoutRefs();
 375       new((*fVertices)[nvertices++]) AliAODVertex(*tmp);
 376       delete tmp;
 377     }
 378   }
 379
 380   fDimuons->Clear("C");
 381
 382   // as there might be a track cut (different from the one of the main filtering),
 383   // we must recreate the dimuon completely from scratch to be 100% safe...
 384
 385   Int_t ndimuons(0);
 386
 387   for ( Int_t i = 0; i < ntracks; ++i )
 388   {
 389     for ( Int_t j = i+1; j < ntracks; ++j )
 390     {
 391       new((*fDimuons)[ndimuons++]) AliAODDimuon(fTracks->At(i),fTracks->At(j));
 392     }
 393   }
 394
 395   AliDebug(1,Form("input mu tracks=%d tracks=%d vertices=%d ndimuons=%d",
 396                   input,fTracks->GetEntries(),fVertices->GetEntries(),fDimuons->GetEntries()));
 397
 398   // Finally, deal with MC information, if needed
 399
 400   if ( fMCMode > 0 )
 401   {
 402     FilterMC(source);
 403   }
 404 }
 405