PWG/muon/AliAODMuonReplicator.cxx

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