PWG/muon/AliAnalysisTaskESDMuonFilter.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 /// Add the muon tracks to the generic AOD track branch during the
  20 /// filtering of the ESD.
  21 ///
  22 ///
  23 /// Note that we :
  24 ///
  25 ///   - completely disable all the branches that are not required by (most) the muon analyses,
  26 ///     e.g. cascades, v0s, kinks, jets, etc...
  27 ///   - filter severely the tracks (keep only muon tracks) and vertices (keep only primary -including
  28 ///     pile-up - vertices) branches
  29 ///
  30 /// \see AliAODMuonReplicator
  31 ///
  32 /// (see AddFilteredAOD method)
  33 ///
  34
  35 #include "AliAnalysisTaskESDMuonFilter.h"
  36
  37 #include "AliAnalysisFilter.h"
  38 #include "AliAnalysisManager.h"
  39 #include "AliAnalysisNonMuonTrackCuts.h"
  40 #include "AliAnalysisNonPrimaryVertices.h"
  41 #include "AliAODDimuon.h"
  42 #include "AliAODEvent.h"
  43 #include "AliAODExtension.h"
  44 #include "AliAODHandler.h"
  45 #include "AliAODMCParticle.h"
  46 #include "AliAODMuonReplicator.h"
  47 #include "AliAODVertex.h"
  48 #include "AliCodeTimer.h"
  49 #include "AliESDEvent.h"
  50 #include "AliESDInputHandler.h"
  51 #include "AliESDMuonTrack.h"
  52 #include "AliESDtrack.h"
  53 #include "AliESDMuonGlobalTrack.h"   // AU
  54 #include "AliMFTAnalysisTools.h"     // AU
  55 #include "AliESDVertex.h"
  56 #include "AliLog.h"
  57 #include "AliMCEvent.h"
  58 #include "AliMCEventHandler.h"
  59 #include "AliMultiplicity.h"
  60 #include <TChain.h>
  61 #include <TFile.h>
  62 #include <TParticle.h>
  63
  64 using std::cout;
  65 using std::endl;
  66 ClassImp(AliAnalysisTaskESDMuonFilter)
  67
  68 AliAnalysisTaskESDMuonFilter::AliAnalysisTaskESDMuonFilter(Bool_t onlyMuon, Bool_t keepAllEvents, Int_t mcMode, Bool_t withSPDtracklets):
  69   AliAnalysisTaskSE(),
  70   fTrackFilter(0x0),
  71   fEnableMuonAOD(kFALSE),
  72   fEnableDimuonAOD(kFALSE),
  73   fOnlyMuon(onlyMuon),
  74   fKeepAllEvents(keepAllEvents),
  75   fMCMode(mcMode),
  76   fWithSPDTracklets(withSPDtracklets)
  77 {
  78   /// Default constructor
  79 }
  80
  81 AliAnalysisTaskESDMuonFilter::AliAnalysisTaskESDMuonFilter(const char* name, Bool_t onlyMuon, Bool_t keepAllEvents, Int_t mcMode, Bool_t withSPDtracklets):
  82   AliAnalysisTaskSE(name),
  83   fTrackFilter(0x0),
  84   fEnableMuonAOD(kFALSE),
  85   fEnableDimuonAOD(kFALSE),
  86   fOnlyMuon(onlyMuon),
  87   fKeepAllEvents(keepAllEvents),
  88   fMCMode(mcMode),
  89   fWithSPDTracklets(withSPDtracklets)
  90 {
  91   /// Constructor
  92 }
  93
  94 //______________________________________________________________________________
  95 void AliAnalysisTaskESDMuonFilter::UserCreateOutputObjects()
  96 {
  97   /// Create the output container
  98   if (fTrackFilter) OutputTree()->GetUserInfo()->Add(fTrackFilter);
  99 }
 100
 101 //______________________________________________________________________________
 102 void AliAnalysisTaskESDMuonFilter::PrintTask(Option_t *option, Int_t indent) const
 103 {
 104   /// Specify how we are configured
 105
 106   AliAnalysisTaskSE::PrintTask(option,indent);
 107
 108   TString spaces(' ',indent+3);
 109
 110   if ( fOnlyMuon )
 111   {
 112     cout << spaces.Data() << "Keep only muon information " << endl;
 113   }
 114   else
 115   {
 116     cout << spaces.Data() << "Keep all information from standard AOD" << endl;
 117   }
 118
 119   if ( fKeepAllEvents )
 120   {
 121     cout << spaces.Data() << "Keep all events, regardless of number of muons" << endl;
 122   }
 123   else
 124   {
 125     cout << spaces.Data() << "Keep only events with at least one muon" << endl;
 126   }
 127
 128   if ( fMCMode > 0 )
 129   {
 130     cout << spaces.Data() << "Assuming work on MC data (i.e. will transmit MC branches)" << endl;
 131     if ( fMCMode == 1 )
 132     {
 133       cout << spaces.Data() << "  (will write MC information irrespective of whether or not we have reconstructed muons in the event)" << endl;
 134     }
 135     else
 136     {
 137       cout << spaces.Data() << "  (will write MC information only if we have reconstructed muons in the event)" << endl;
 138     }
 139   }
 140
 141   if ( fWithSPDTracklets )
 142   {
 143     cout << spaces.Data() << "Will also keep SPD tracklets" << endl;
 144   }
 145 }
 146
 147 //______________________________________________________________________________
 148 void AliAnalysisTaskESDMuonFilter::AddFilteredAOD(const char* aodfilename, const char* title)
 149 {
 150   /// Add an output filtered and replicated aod
 151
 152   AliAODHandler *aodH = (AliAODHandler*)((AliAnalysisManager::GetAnalysisManager())->GetOutputEventHandler());
 153   if (!aodH) Fatal("UserCreateOutputObjects", "No AOD handler");
 154
 155   AliAODExtension* ext = aodH->AddFilteredAOD(aodfilename,title);
 156
 157   if (!ext) return;
 158
 159   if ( fOnlyMuon )
 160   {
 161
 162     AliAODMuonReplicator* murep = new AliAODMuonReplicator("MuonReplicator",
 163                                                            "remove non muon tracks and non primary or pileup vertices",
 164                                                            new AliAnalysisNonMuonTrackCuts,
 165                                                            new AliAnalysisNonPrimaryVertices,
 166                                                            fMCMode,
 167                                                            kFALSE,
 168                                                            fWithSPDTracklets);
 169
 170     ext->DropUnspecifiedBranches(); // all branches not part of a FilterBranch call (below) will be dropped
 171
 172     ext->FilterBranch("tracks",murep);
 173     ext->FilterBranch("vertices",murep);
 174     ext->FilterBranch("dimuons",murep);
 175     ext->FilterBranch("AliAODVZERO",murep);
 176     ext->FilterBranch("AliAODTZERO",murep);
 177     ext->FilterBranch("AliAODZDC",murep);
 178
 179     if ( fWithSPDTracklets )
 180     {
 181       ext->FilterBranch("tracklets",murep);
 182     }
 183
 184     if ( fMCMode > 0 )
 185     {
 186       // MC branches will be copied (if present), as they are, but only
 187       // for events with at least one muon.
 188       // For events w/o muon, mcparticles array will be empty and mcheader will be dummy
 189       // (e.g. strlen(GetGeneratorName())==0)
 190
 191       ext->FilterBranch("mcparticles",murep);
 192       ext->FilterBranch("mcHeader",murep);
 193     }
 194   }
 195 }
 196
 197 //______________________________________________________________________________
 198 void AliAnalysisTaskESDMuonFilter::Init()
 199 {
 200   /// Initialization
 201   if(fEnableMuonAOD) AddFilteredAOD("AliAOD.Muons.root", "MuonEvents");
 202   if(fEnableDimuonAOD) AddFilteredAOD("AliAOD.Dimuons.root", "DimuonEvents");
 203 }
 204
 205
 206 //______________________________________________________________________________
 207 void AliAnalysisTaskESDMuonFilter::UserExec(Option_t */*option*/)
 208 {
 209   /// Execute analysis for current event
 210
 211   Long64_t ientry = Entry();
 212   if(fDebug)printf("Muon Filter: Analysing event # %5d\n", (Int_t) ientry);
 213
 214   ConvertESDtoAOD();
 215 }
 216
 217 //______________________________________________________________________________
 218 void AliAnalysisTaskESDMuonFilter::ConvertESDtoAOD()
 219 {
 220   /// ESD Muon Filter analysis task executed for each event
 221
 222   AliCodeTimerAuto("",0);
 223
 224   AliESDEvent* esd = dynamic_cast<AliESDEvent*>(InputEvent());
 225   if (!esd)
 226   {
 227     AliError("Could not get input ESD event");
 228     return;
 229   }
 230
 231   AliMCEventHandler *mcH = static_cast<AliMCEventHandler*>((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
 232
 233   // Define arrays for muons
 234   Double_t pos[3];
 235   Double_t p[3];
 236   //  Double_t pid[10];
 237
 238   // has to be changed once the muon pid is provided by the ESD
 239   //  for (Int_t i = 0; i < 10; pid[i++] = 0.) {}
 240   //  pid[AliAODTrack::kMuon]=1.;
 241
 242   AliAODHeader* header = AODEvent()->GetHeader();
 243   AliAODTrack *aodTrack = 0x0;
 244   AliESDMuonTrack *esdMuTrack = 0x0;
 245
 246   // Access to the AOD container of tracks
 247   TClonesArray &tracks = *(AODEvent()->GetTracks());
 248   Int_t jTracks = tracks.GetEntriesFast();
 249
 250   // Read primary vertex from AOD event
 251   AliAODVertex *primary = AODEvent()->GetPrimaryVertex();
 252   if (fDebug && primary) primary->Print();
 253
 254   // ----------------- MUON TRACKS -----------------------------------------------------
 255
 256   // Loop on muon tracks to fill the AOD track branch
 257   Int_t nMuTracks = esd->GetNumberOfMuonTracks();
 258
 259   for (Int_t iTrack=0; iTrack<nMuTracks; ++iTrack) esd->GetMuonTrack(iTrack)->SetESDEvent(esd);
 260
 261   // Update number of positive and negative tracks from AOD event (M.G.)
 262   Int_t nTracks    = header->GetRefMultiplicity();
 263   Int_t nPosTracks = header->GetRefMultiplicityPos();
 264   Int_t nNegTracks = header->GetRefMultiplicityNeg();
 265
 266   // Access to the AOD container of dimuons
 267   TClonesArray &dimuons = *(AODEvent()->GetDimuons());
 268
 269   Int_t nMuons=0;
 270   Int_t nDimuons=0;
 271   Int_t jDimuons=0;
 272   Int_t nMuonTrack[100];
 273   UChar_t itsClusMap(0);
 274
 275   for(int imuon=0;imuon<100;imuon++) nMuonTrack[imuon]=0;
 276
 277   for (Int_t nMuTrack = 0; nMuTrack < nMuTracks; ++nMuTrack)
 278   {
 279     esdMuTrack = esd->GetMuonTrack(nMuTrack);
 280
 281     if (!esdMuTrack->ContainTrackerData()) continue;
 282
 283     UInt_t selectInfo(0);
 284
 285     // Track selection
 286     if (fTrackFilter) {
 287         selectInfo = fTrackFilter->IsSelected(esdMuTrack);
 288         if (!selectInfo) {
 289           continue;
 290         }
 291     }
 292
 293     p[0] = esdMuTrack->Px();
 294     p[1] = esdMuTrack->Py();
 295     p[2] = esdMuTrack->Pz();
 296
 297     pos[0] = esdMuTrack->GetNonBendingCoor();
 298     pos[1] = esdMuTrack->GetBendingCoor();
 299     pos[2] = esdMuTrack->GetZ();
 300
 301     if (mcH) mcH->SelectParticle(esdMuTrack->GetLabel()); // to insure that particle's ancestors will be in output MC branches
 302
 303     aodTrack = new(tracks[jTracks++]) AliAODTrack(esdMuTrack->GetUniqueID(), // ID
 304                                                   esdMuTrack->GetLabel(), // label
 305                                                   p, // momentum
 306                                                   kTRUE, // cartesian coordinate system
 307                                                   pos, // position
 308                                                   kFALSE, // isDCA
 309                                                   0x0, // covariance matrix
 310                                                   esdMuTrack->Charge(), // charge
 311                                                   itsClusMap, // ITSClusterMap
 312                                                   //pid, // pid
 313                                                   primary, // primary vertex
 314                                                   kFALSE, // used for vertex fit?
 315                                                   kFALSE, // used for primary vertex fit?
 316                                                   AliAODTrack::kPrimary,// track type
 317                                                   selectInfo);
 318
 319     aodTrack->SetPIDForTracking(AliPID::kMuon);
 320     aodTrack->SetXYAtDCA(esdMuTrack->GetNonBendingCoorAtDCA(), esdMuTrack->GetBendingCoorAtDCA());
 321     aodTrack->SetPxPyPzAtDCA(esdMuTrack->PxAtDCA(), esdMuTrack->PyAtDCA(), esdMuTrack->PzAtDCA());
 322     aodTrack->SetRAtAbsorberEnd(esdMuTrack->GetRAtAbsorberEnd());
 323     aodTrack->ConvertAliPIDtoAODPID();
 324     aodTrack->SetChi2perNDF(esdMuTrack->GetChi2() / (2.*esdMuTrack->GetNHit() - 5.));
 325     aodTrack->SetChi2MatchTrigger(esdMuTrack->GetChi2MatchTrigger());
 326     aodTrack->SetHitsPatternInTrigCh(esdMuTrack->GetHitsPatternInTrigCh());
 327     UInt_t pattern = esdMuTrack->GetHitsPatternInTrigCh();
 328     AliESDMuonTrack::AddEffInfo(pattern, 0, esdMuTrack->LoCircuit(), (AliESDMuonTrack::EAliTriggerChPatternFlag)0);
 329     aodTrack->SetMUONtrigHitsMapTrg(pattern);
 330     aodTrack->SetMUONtrigHitsMapTrk(esdMuTrack->GetHitsPatternInTrigChTrk());
 331     aodTrack->SetMuonClusterMap(esdMuTrack->GetMuonClusterMap());
 332     aodTrack->SetMatchTrigger(esdMuTrack->GetMatchTrigger());
 333     aodTrack->Connected(esdMuTrack->IsConnected());
 334     primary->AddDaughter(aodTrack);
 335
 336     ++nTracks;
 337     if (esdMuTrack->Charge() > 0) nPosTracks++;
 338     else nNegTracks++;
 339
 340     nMuonTrack[nMuons]= jTracks-1;
 341     ++nMuons;
 342   }
 343
 344   if(nMuons>=2)
 345   {
 346     for(int i=0;i<nMuons;i++){
 347       Int_t index0 = nMuonTrack[i];
 348       for(int j=i+1;j<nMuons;j++){
 349         Int_t index1 = nMuonTrack[j];
 350         new(dimuons[jDimuons++]) AliAODDimuon(tracks.At(index0),tracks.At(index1));
 351         ++nDimuons;
 352         if (fDebug > 1){
 353           AliAODDimuon *dimuon0 = (AliAODDimuon*)dimuons.At(jDimuons-1);
 354           printf("Dimuon: mass = %f, px=%f, py=%f, pz=%f\n",dimuon0->M(),dimuon0->Px(),dimuon0->Py(),dimuon0->Pz());
 355           AliAODTrack  *mu0 = (AliAODTrack*) dimuon0->GetMu(0);
 356           AliAODTrack  *mu1 = (AliAODTrack*) dimuon0->GetMu(1);
 357           printf("Muon0 px=%f py=%f pz=%f\n",mu0->Px(),mu0->Py(),mu0->Pz());
 358           printf("Muon1 px=%f py=%f pz=%f\n",mu1->Px(),mu1->Py(),mu1->Pz());
 359         }
 360       }
 361     }
 362   }
 363
 364
 365   header->SetRefMultiplicity(nTracks);
 366   header->SetRefMultiplicityPos(nPosTracks);
 367   header->SetRefMultiplicityNeg(nNegTracks);
 368   header->SetNumberOfMuons(nMuons);
 369   header->SetNumberOfDimuons(nDimuons);
 370
 371   // ----------------- MFT + MUON TRACKS -----------------------------------------------------
 372
 373   AliESDMuonGlobalTrack *esdMuGlobalTrack = 0x0;
 374
 375   // Loop on muon global tracks to fill the AOD track branch.
 376   // We won't update the number of total, pos and neg tracks in the event
 377
 378   Int_t nMuGlobalTracks = esd->GetNumberOfMuonGlobalTracks();
 379
 380   for (Int_t iTrack=0; iTrack<nMuGlobalTracks; ++iTrack) esd->GetMuonGlobalTrack(iTrack)->SetESDEvent(esd);
 381
 382   Int_t nGlobalMuons=0;
 383   Int_t nGlobalDimuons=0;
 384   Int_t jGlobalDimuons=0;
 385   Int_t nMuonGlobalTrack[100];
 386   itsClusMap = 0;
 387
 388   for (Int_t iMuon=0; iMuon<100; iMuon++) nMuonGlobalTrack[iMuon]=0;    // position of the i-th muon track in the tracks array of the AOD event
 389
 390   for (Int_t nMuTrack=0; nMuTrack<nMuGlobalTracks; ++nMuTrack) {
 391
 392     esdMuGlobalTrack = esd->GetMuonGlobalTrack(nMuTrack);
 393
 394     if (!esdMuGlobalTrack->ContainTrackerData()) continue;
 395
 396     UInt_t selectInfo(0);
 397
 398     // Track selection
 399     if (fTrackFilter) {
 400       selectInfo = fTrackFilter->IsSelected(esdMuGlobalTrack);
 401       if (!selectInfo) {
 402         continue;
 403       }
 404     }
 405
 406     p[0] = esdMuGlobalTrack->Px();
 407     p[1] = esdMuGlobalTrack->Py();
 408     p[2] = esdMuGlobalTrack->Pz();
 409
 410     esdMuGlobalTrack -> GetFirstTrackingPoint(pos);
 411
 412     if (mcH) mcH->SelectParticle(esdMuGlobalTrack->GetLabel()); // to insure that particle's ancestors will be in output MC branches
 413
 414     Double_t covTr[21] = {0};
 415     AliMFTAnalysisTools::ConvertCovMatrixMUON2AOD(esdMuGlobalTrack->GetCovariances(), covTr);
 416
 417     aodTrack = new(tracks[jTracks++]) AliAODTrack(esdMuGlobalTrack->GetUniqueID(), // ID
 418                                                   esdMuGlobalTrack->GetLabel(),    // label
 419                                                   p,                               // momentum
 420                                                   kTRUE,                           // cartesian coordinate system
 421                                                   pos,                             // position
 422                                                   kFALSE,                          // isDCA
 423                                                   covTr,                           // covariance matrix
 424                                                   esdMuGlobalTrack->Charge(),      // charge
 425                                                   itsClusMap,                      // ITSClusterMap
 426                                                   primary,                         // origin vertex
 427                                                   kFALSE,                          // used for vertex fit?
 428                                                   kFALSE,                          // used for primary vertex fit?
 429                                                   AliAODTrack::kPrimary,           // track type
 430                                                   selectInfo);
 431
 432
 433     aodTrack->SetPIDForTracking(AliPID::kMuon);
 434
 435     Double_t xyAtVertex[2] = {0};
 436     esdMuGlobalTrack -> GetXYAtVertex(xyAtVertex);
 437
 438     aodTrack->SetIsMuonGlobalTrack(kTRUE);
 439
 440     aodTrack->SetMFTClusterPattern(esdMuGlobalTrack->GetMFTClusterPattern());
 441     aodTrack->SetXYAtDCA(xyAtVertex[0], xyAtVertex[1]);
 442     aodTrack->SetPxPyPzAtDCA(p[0], p[1], p[2]);
 443     aodTrack->SetRAtAbsorberEnd(esdMuGlobalTrack->GetRAtAbsorberEnd());
 444     aodTrack->ConvertAliPIDtoAODPID();
 445     aodTrack->SetChi2perNDF(esdMuGlobalTrack->GetChi2OverNdf());
 446     aodTrack->SetChi2MatchTrigger(esdMuGlobalTrack->GetChi2MatchTrigger());
 447     aodTrack->SetHitsPatternInTrigCh(esdMuGlobalTrack->GetHitsPatternInTrigCh());
 448     UInt_t pattern = esdMuGlobalTrack->GetHitsPatternInTrigCh();
 449     AliESDMuonTrack::AddEffInfo(pattern, 0, esdMuGlobalTrack->GetLoCircuit(), (AliESDMuonTrack::EAliTriggerChPatternFlag)0);
 450     aodTrack->SetMUONtrigHitsMapTrg(pattern);
 451     aodTrack->SetMUONtrigHitsMapTrk(esdMuGlobalTrack->GetHitsPatternInTrigChTrk());
 452     aodTrack->SetMuonClusterMap(esdMuGlobalTrack->GetMuonClusterMap());
 453     aodTrack->SetMatchTrigger(esdMuGlobalTrack->GetMatchTrigger());
 454     aodTrack->Connected(esdMuGlobalTrack->IsConnected());
 455
 456     printf("Added Muon Global Track %d of %d\n",nMuTrack,nMuGlobalTracks);
 457     aodTrack->Print();     // to be removed!!!
 458
 459     primary->AddDaughter(aodTrack);
 460
 461     nMuonGlobalTrack[nGlobalMuons] = jTracks-1;
 462     ++nGlobalMuons;
 463   }
 464
 465   if (nGlobalMuons >= 2) {
 466     for (Int_t i=0; i<nGlobalMuons; i++) {
 467       Int_t index0 = nMuonGlobalTrack[i];
 468       for (Int_t j=i+1; j<nGlobalMuons; j++) {
 469         Int_t index1 = nMuonGlobalTrack[j];
 470         new (dimuons[jGlobalDimuons++]) AliAODDimuon(tracks.At(index0), tracks.At(index1));
 471         ++nDimuons;
 472         if (fDebug > 1) {
 473           AliAODDimuon *dimuon0 = (AliAODDimuon*)dimuons.At(jGlobalDimuons-1);
 474           printf("Dimuon: mass = %f, px=%f, py=%f, pz=%f\n",dimuon0->M(),dimuon0->Px(),dimuon0->Py(),dimuon0->Pz());
 475           AliAODTrack  *mu0 = (AliAODTrack*) dimuon0->GetMu(0);
 476           AliAODTrack  *mu1 = (AliAODTrack*) dimuon0->GetMu(1);
 477           printf("Muon0 px=%f py=%f pz=%f\n",mu0->Px(),mu0->Py(),mu0->Pz());
 478           printf("Muon1 px=%f py=%f pz=%f\n",mu1->Px(),mu1->Py(),mu1->Pz());
 479         }
 480       }
 481     }
 482   }
 483
 484   header->SetNumberOfGlobalMuons(nGlobalMuons);
 485   header->SetNumberOfGlobalDimuons(nGlobalDimuons);
 486
 487   // -----------------------------------------------------------------------
 488
 489   AliAODHandler* handler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 490
 491   if ( handler && fEnableMuonAOD && ( (nMuons>0) || fKeepAllEvents ) )
 492   {
 493     AliAODExtension *extMuons = handler->GetFilteredAOD("AliAOD.Muons.root");
 494     if ( extMuons ) extMuons->SelectEvent();
 495   }
 496
 497   if ( handler && fEnableDimuonAOD && ( (nMuons>1) || fKeepAllEvents )  )
 498   {
 499     AliAODExtension *extDimuons = handler->GetFilteredAOD("AliAOD.Dimuons.root");
 500     if ( extDimuons ) extDimuons->SelectEvent();
 501   }
 502
 503 }