PWGUD/UPC/AliAnalysisTaskUpcFilterSemiforward.cxx

   1
   2 //_____________________________________________________________________________
   3 //    Class for semiforward UPC filter
   4 //    Author: Jaroslav Adam
   5 //
   6 //    Fill structure of AliUPCEvent
   7 //_____________________________________________________________________________
   8
   9 // c++ headers
  10 #include <iostream>
  11 #include <string.h>
  12
  13 // root headers
  14 #include "TList.h"
  15 #include "TH1I.h"
  16 #include "TTree.h"
  17 #include "TClonesArray.h"
  18 #include "TParticle.h"
  19 #include "TObjString.h"
  20 #include "TFile.h"
  21 #include "TArrayF.h"
  22
  23 // aliroot headers
  24 #include "AliAnalysisManager.h"
  25 #include "AliInputEventHandler.h"
  26 #include "AliESDEvent.h"
  27 #include "AliAODEvent.h"
  28 #include "AliMCEvent.h"
  29 #include "AliAODVZERO.h"
  30 #include "AliAODZDC.h"
  31 #include "AliAODVertex.h"
  32 #include "AliESDVZERO.h"
  33 #include "AliESDZDC.h"
  34 #include "AliESDVertex.h"
  35 #include "AliAODTrack.h"
  36 #include "AliMultiplicity.h"
  37 #include "AliESDtrack.h"
  38 #include "AliESDMuonTrack.h"
  39 #include "AliAODMCParticle.h"
  40 #include "AliMCParticle.h"
  41 #include "AliTriggerAnalysis.h"
  42 #include "AliMuonTrackCuts.h"
  43 #include "AliESDtrackCuts.h"
  44 #include "AliAODPid.h"
  45 #include "AliGenEventHeader.h"
  46 #include "AliAODMCHeader.h"
  47
  48 // my headers
  49 #include "AliAnalysisTaskUpcFilterSemiforward.h"
  50 #include "AliUPCTrack.h"
  51 #include "AliUPCMuonTrack.h"
  52 #include "AliUPCEvent.h"
  53
  54 ClassImp(AliAnalysisTaskUpcFilterSemiforward);
  55
  56 // task for upc semiforward filter
  57 // jaroslav.adam@cern.ch
  58
  59 //_____________________________________________________________________________
  60 AliAnalysisTaskUpcFilterSemiforward::AliAnalysisTaskUpcFilterSemiforward(const char *name)
  61  :AliAnalysisTaskSE(name),
  62   fIsESD(0), fIsMC(0), fMuonCuts(0x0), fTriggerAna(0x0), fCutsList(0x0),
  63   fHistList(0x0), fCounter(0x0), fTriggerCounter(0x0),
  64   fUPCEvent(0x0), fUPCTree(0x0)
  65 {
  66
  67   // Constructor
  68
  69   DefineOutput(1, TTree::Class());
  70   DefineOutput(2, TList::Class());
  71
  72 }//AliAnalysisTaskUpcFilterSemiforward
  73
  74 //_____________________________________________________________________________
  75 AliAnalysisTaskUpcFilterSemiforward::~AliAnalysisTaskUpcFilterSemiforward()
  76 {
  77   // destructor
  78
  79   if(fHistList) {delete fHistList; fHistList = 0x0;}
  80   if(fCounter) {delete fCounter; fCounter = 0x0;}
  81   if(fTriggerCounter) {delete fTriggerCounter; fTriggerCounter = 0x0;}
  82   if(fUPCEvent) {delete fUPCEvent; fUPCEvent = 0x0;}
  83   if(fUPCTree) {delete fUPCTree; fUPCTree = 0x0;}
  84   if(fMuonCuts) {delete fMuonCuts; fMuonCuts = 0x0;}
  85   if(fTriggerAna) {delete fTriggerAna; fTriggerAna = 0x0;}
  86   if(fCutsList) {delete[] fCutsList; fCutsList=0x0;}
  87
  88 }//~AliAnalysisTaskUpcFilterSemiforward
  89
  90 //_____________________________________________________________________________
  91 void AliAnalysisTaskUpcFilterSemiforward::UserCreateOutputObjects()
  92 {
  93   //muon track cuts
  94   fMuonCuts = new AliMuonTrackCuts("StdMuonCuts", "StdMuonCuts");
  95   fMuonCuts->SetFilterMask ( AliMuonTrackCuts::kMuPdca );
  96   fMuonCuts->Print("mask");
  97   fMuonCuts->SetAllowDefaultParams(kTRUE);
  98   fMuonCuts->SetPassName("muon_pass2");
  99
 100   //trigger analysis for SPD FO fired chips in MC
 101   fTriggerAna = new AliTriggerAnalysis();
 102   fTriggerAna->SetAnalyzeMC( fIsMC );
 103
 104   //ESD track cuts
 105   if( fIsESD ) {
 106     fCutsList = new AliESDtrackCuts*[32];
 107     for(UInt_t i=0; i<32; i++) fCutsList[i] = 0x0;
 108
 109     // bits 0 - 10 set in $ALICE_ROOT/ANALYSIS/ESDfilter/macros/AddTaskESDFilter.C
 110
 111     // Cuts on primary tracks
 112     AliESDtrackCuts* esdTrackCutsL = AliESDtrackCuts::GetStandardTPCOnlyTrackCuts();
 113     fCutsList[0] = esdTrackCutsL;
 114
 115     // standard cuts with very loose DCA
 116     AliESDtrackCuts* esdTrackCutsH = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011(kFALSE);
 117     esdTrackCutsH->SetMaxDCAToVertexXY(2.4);
 118     esdTrackCutsH->SetMaxDCAToVertexZ(3.2);
 119     esdTrackCutsH->SetDCAToVertex2D(kTRUE);
 120     fCutsList[4] = esdTrackCutsH;
 121
 122     // standard cuts with tight DCA cut
 123     AliESDtrackCuts *esdTrackCutsH2 = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011();
 124     fCutsList[5] = esdTrackCutsH2;
 125
 126     // bits 11 - 31 free
 127
 128     // selection by Daniele, UPC meeting, 15 July 2014
 129     // https://indico.cern.ch/event/330483/contribution/3/material/slides/0.pdf
 130
 131     AliESDtrackCuts *cuts11 = new AliESDtrackCuts();
 132     cuts11->SetMinNClustersTPC(70);
 133     cuts11->SetRequireTPCRefit(kTRUE);
 134     cuts11->SetRequireITSRefit(kTRUE);
 135     cuts11->SetMaxChi2PerClusterTPC(4);
 136     cuts11->SetAcceptKinkDaughters(kFALSE);
 137     cuts11->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
 138     cuts11->SetMaxDCAToVertexXYPtDep("0.0105+0.0350/pt^1.01");
 139     cuts11->SetMaxDCAToVertexZ(2);
 140     cuts11->SetDCAToVertex2D(kFALSE);
 141     cuts11->SetRequireSigmaToVertex(kFALSE);
 142     cuts11->SetMaxChi2PerClusterITS(5.);
 143     fCutsList[11] = cuts11;
 144
 145     AliESDtrackCuts *cuts12 = new AliESDtrackCuts();
 146     cuts12->SetMinNClustersTPC(70);
 147     cuts12->SetRequireTPCRefit(kTRUE);
 148     cuts12->SetRequireITSRefit(kTRUE);
 149     cuts12->SetMaxChi2PerClusterTPC(4);
 150     cuts12->SetAcceptKinkDaughters(kFALSE);
 151     cuts12->SetClusterRequirementITS(AliESDtrackCuts::kSPD,AliESDtrackCuts::kAny);
 152     cuts12->SetMaxDCAToVertexXYPtDep("0.0105+0.0350/pt^1.01");
 153     cuts12->SetMaxDCAToVertexZ(2);
 154     cuts12->SetDCAToVertex2D(kFALSE);
 155     cuts12->SetRequireSigmaToVertex(kFALSE);
 156     fCutsList[12] = cuts12;
 157
 158     // setting in Evgeny's trees
 159     AliESDtrackCuts *cuts13 = AliESDtrackCuts::GetStandardITSTPCTrackCuts2011();
 160     cuts13->SetMaxChi2TPCConstrainedGlobal(1e+10);
 161     fCutsList[13] = cuts13;
 162
 163   }
 164
 165   //output histograms
 166   fHistList = new TList();
 167   fHistList->SetOwner();
 168   fCounter = new TH1I("fCounter", "fCounter", 30, 1, 31);
 169   fHistList->Add(fCounter);
 170   fTriggerCounter = new TH2I("fTriggerCounter", "fTriggerCounter", 3000, 195000, 198000, NTRG+1, 0, NTRG+1);
 171   fHistList->Add(fTriggerCounter);
 172
 173   //output tree
 174   fUPCEvent  = new AliUPCEvent();
 175   fUPCEvent->SetIsMC( fIsMC );
 176   fUPCEvent->SetIsESD( fIsESD );
 177
 178   TDirectory *pwd = gDirectory;
 179   OpenFile(1);
 180   fUPCTree = new TTree("fUPCTree", "fUPCTree");
 181   pwd->cd();
 182   fUPCTree->Branch("fUPCEvent", &fUPCEvent);
 183
 184   PostData(1, fUPCTree);
 185   PostData(2, fHistList);
 186
 187 }//UserCreateOutputObjects
 188
 189 //_____________________________________________________________________________
 190 void AliAnalysisTaskUpcFilterSemiforward::NotifyRun()
 191 {
 192
 193   //cout<<"AliAnalysisTaskUpcFilterSemiforward::NotifyRun()"<<endl;
 194
 195   fMuonCuts->SetRun(fInputHandler); // use input handler from AliAnalysisTaskSE
 196
 197   if( fIsESD ) { ((AliESDEvent*) InputEvent())->InitMagneticField(); }
 198
 199 }//NotifyRun
 200
 201 //_____________________________________________________________________________
 202 void AliAnalysisTaskUpcFilterSemiforward::UserExec(Option_t *)
 203 {
 204
 205   //cout<<"#################### Next event ##################"<<endl;
 206
 207   // input event
 208   AliVEvent *vEvent = InputEvent();
 209   if(!vEvent) return;
 210
 211   fUPCEvent->ClearEvent();
 212
 213   fCounter->Fill( 1 ); // 1 = analyzed events
 214
 215   // trigger
 216   TString trigger = vEvent->GetFiredTriggerClasses();
 217   Bool_t trgClasses[NTRG]; // array of fired trigger classes
 218   for(Int_t itrg=0; itrg<NTRG; itrg++) trgClasses[itrg] = kFALSE;
 219
 220   trgClasses[1] = trigger.Contains("CMUP6-B"); // p-Pb FW
 221   trgClasses[2] = trigger.Contains("CMUP3-B"); // Pb-p FW
 222   trgClasses[3] = trigger.Contains("CMUP8-B"); // Pb-p FW
 223
 224   trgClasses[4] = trigger.Contains("CMUP7-B"); // p-Pb SFW
 225   trgClasses[5] = trigger.Contains("CMUP5-B"); // Pb-p SFW
 226   trgClasses[6] = trigger.Contains("CMUP9-B"); // Pb-p SFW
 227
 228   trgClasses[7] = trigger.Contains("CCUP7-B"); // CEN
 229
 230   trgClasses[8] = trigger.Contains("CMUP7-ACE");
 231   trgClasses[9] = trigger.Contains("CMUP5-ACE");
 232   trgClasses[10]= trigger.Contains("CMUP9-ACE");
 233
 234   Bool_t isTrg = kFALSE;
 235   for(Int_t itrg=1; itrg<NTRG; itrg++) {
 236     if(!trgClasses[itrg]) continue;
 237     //trigger at itrg is fired
 238     fUPCEvent->SetTriggerClass( itrg , kTRUE );
 239     fTriggerCounter->Fill( vEvent->GetRunNumber() , itrg );
 240     isTrg = kTRUE;
 241   }
 242   if(!isTrg && !fIsMC) {PostData(2, fHistList); return;}
 243   //event passed the trigger
 244
 245   fCounter->Fill( 2 ); // 2 = events after trigger (ESD and AOD)
 246
 247   // ESD / AOD specific tasks: MC, SPD FO, L0 inputs, SPD tracklets, tracks, ZDC tdc
 248   Bool_t stat;
 249   if( fIsESD ) stat = RunESD();
 250   if(!fIsESD ) stat = RunAOD();
 251   if( !stat && !fIsMC ) {PostData(2, fHistList); return;}
 252   //event passed ESD / AOD specific selection
 253
 254   fCounter->Fill( 3 ); // 3 = events after ESD / AOD specific part
 255
 256   // input data
 257   const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
 258   // reconstruction pass: -1 = unknown, 1 = pass1, 2 = pass2, counter indices: 9 (unknown), 11 (pass1), 12 (pass2)
 259   fUPCEvent->SetRecoPass( -1 );
 260   if( strstr(filnam,"/pass1/") ) {fUPCEvent->SetRecoPass( 1 ); fCounter->Fill( 11 );}
 261   if( strstr(filnam,"/pass2/") ) {fUPCEvent->SetRecoPass( 2 ); fCounter->Fill( 12 );}
 262   if( fUPCEvent->GetRecoPass() < 0 ) fCounter->Fill( 9 );
 263   fUPCEvent->SetInputFileName( filnam );
 264   fUPCEvent->SetEventNumber( ((TTree*) GetInputData(0))->GetTree()->GetReadEntry() );
 265   fUPCEvent->SetRunNumber( vEvent->GetRunNumber() );
 266
 267   //VZERO
 268   AliVVZERO *dataVZERO = vEvent->GetVZEROData();
 269   if(!dataVZERO) {PostData(2, fHistList); return;}
 270
 271   fUPCEvent->SetV0ADecision( dataVZERO->GetV0ADecision() );
 272   fUPCEvent->SetV0CDecision( dataVZERO->GetV0CDecision() );
 273   for (UInt_t iv=0; iv<32; iv++) {
 274     if( dataVZERO->BBTriggerV0C((Int_t)iv) ) fUPCEvent->SetBBtriggerV0Cmask(iv);
 275     if( dataVZERO->GetBBFlag((Int_t)iv) ) fUPCEvent->SetBBFlagV0Cmask(iv);
 276   }
 277
 278   //ZDC
 279   AliVZDC *dataZDC = vEvent->GetZDCData();
 280   if(!dataZDC) {PostData(2, fHistList); return;}
 281
 282   //energy in ZDC
 283   Double_t eZnc=0., eZpc=0., eZna=0., eZpa=0.;
 284   for(Int_t i=0; i<5; i++) {
 285     eZnc += dataZDC->GetZNCTowerEnergy()[i];
 286     eZpc += dataZDC->GetZPCTowerEnergy()[i];
 287     eZna += dataZDC->GetZNATowerEnergy()[i];
 288     eZpa += dataZDC->GetZPATowerEnergy()[i];
 289   }
 290   fUPCEvent->SetZNCEnergy( eZnc );
 291   fUPCEvent->SetZPCEnergy( eZpc );
 292   fUPCEvent->SetZNAEnergy( eZna );
 293   fUPCEvent->SetZPAEnergy( eZpa );
 294
 295   //default primary vertex
 296   const AliVVertex *vtx = vEvent->GetPrimaryVertex();
 297   if(!vtx) {PostData(2, fHistList); return;}
 298   Double_t pos[3]; vtx->GetXYZ(pos);
 299   Double_t cov[6]; vtx->GetCovarianceMatrix(cov);
 300   fUPCEvent->SetPrimaryVertex(pos, vtx->GetChi2perNDF(), cov, vtx->GetNContributors());
 301   const char *vtxtitle = vtx->GetTitle();
 302   fUPCEvent->SetPrimaryVertexTitle( vtxtitle );
 303
 304   fCounter->Fill( 4 ); // 4 = events written to the tree (ESD and AOD)
 305
 306   fUPCTree ->Fill();
 307   PostData(1, fUPCTree);
 308   PostData(2, fHistList);
 309
 310 }//UserExec
 311
 312 //_____________________________________________________________________________
 313 Bool_t AliAnalysisTaskUpcFilterSemiforward::RunAOD()
 314 {
 315   //cout<<"#################### AOD event ##################"<<endl;
 316
 317   //input AOD event
 318   AliAODEvent *aodEvent = (AliAODEvent*) InputEvent();
 319   if(!aodEvent) return kFALSE;
 320
 321   fCounter->Fill( 22 ); // 22 = AOD analyzed events
 322
 323   if(fIsMC) RunAODMC( (TClonesArray*) aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName()),
 324                       (AliAODMCHeader*) aodEvent->FindListObject(AliAODMCHeader::StdBranchName())
 325                     );
 326
 327   //nominal interaction point
 328   static AliAODVertex *vtx0 = 0x0;
 329   if( !vtx0 ) {
 330     // make a vertex at coordinates 0, 0, 0
 331     vtx0 = new AliAODVertex();
 332     vtx0->SetX(0.); vtx0->SetY(0.); vtx0->SetZ(0.);
 333   }
 334
 335   //array of tracks for DCAs calculation
 336   //static TClonesArray *trackArray = 0x0;
 337   //if( !trackArray ) trackArray = new TClonesArray("AliAODTrack", 3); // number of DCAs
 338
 339   Double_t pxpypz[3];
 340   UChar_t maskMan;
 341   Double_t xyzDca[2], cov[3];
 342   Float_t b[2], covF[3];
 343   Int_t nmun=0, ncen=0;
 344   // AOD tracks loop
 345   for(Int_t itr=0; itr<aodEvent->GetNumberOfTracks(); itr++) {
 346     AliAODTrack *trk = aodEvent->GetTrack(itr);
 347     if( !trk ) continue;
 348
 349     //muon track
 350     if( trk->IsMuonTrack() ) {
 351
 352       AliUPCMuonTrack *upcMuon = fUPCEvent->AddMuonTrack();
 353       upcMuon->SetPtEtaPhi( trk->Pt(), trk->Eta(), trk->Phi() );
 354       upcMuon->SetCharge( trk->Charge() );
 355       upcMuon->SetMatchTrigger( trk->GetMatchTrigger() );
 356       upcMuon->SetRAtAbsorberEnd( trk->GetRAtAbsorberEnd() );
 357       upcMuon->SetChi2perNDF( trk->Chi2perNDF() );
 358       upcMuon->SetDCA( trk->DCA() );
 359       upcMuon->SetPxDCA( fMuonCuts->IsSelected(trk) );
 360
 361       nmun++;
 362     }
 363     //central track
 364     else {
 365
 366       maskMan = 0;
 367       if( trk->HasPointOnITSLayer(0) )                 { maskMan |= 1 << 0; }
 368       if( trk->HasPointOnITSLayer(1) )                 { maskMan |= 1 << 1; }
 369       if( trk->GetStatus() & AliESDtrack::kITSrefit )  { maskMan |= 1 << 2; }
 370       if( trk->GetStatus() & AliESDtrack::kTPCrefit )  { maskMan |= 1 << 3; }
 371
 372       //selection for at least one point in ITS and TPC refit
 373       //if( !(maskMan & (1 << 0)) && !(maskMan & (1 << 1)) ) continue;
 374       //if( !(maskMan & (1 << 3)) ) continue;
 375
 376       //selection for its refit and tpc refit
 377       //if( !(maskMan & (1 << 2)) || !(maskMan & (1 << 3)) ) continue;
 378
 379       //selection for tpc refit
 380       //if( !(maskMan & (1 << 3)) ) continue;
 381
 382       //selection for at least one point in SPD and ITS refit and TPC refit
 383       //if( !(maskMan & (1 << 0)) && !(maskMan & (1 << 1)) ) continue;
 384       //if( !(maskMan & (1 << 2)) || !(maskMan & (1 << 3)) ) continue;
 385
 386       //selection for at least one point in SPD and at least one TPC cluster
 387       if( !(maskMan & (1 << 0)) && !(maskMan & (1 << 1)) ) continue;
 388       if( trk->GetTPCNcls() == 0 ) continue;
 389
 390       trk->GetPxPyPz(pxpypz);
 391       AliAODPid *apid = trk->GetDetPid();
 392       if(!apid) continue;
 393
 394       AliUPCTrack *upcTrack = fUPCEvent->AddTrack();
 395       upcTrack->SetPxPyPz( pxpypz );
 396       upcTrack->SetMaskMan( maskMan );
 397       upcTrack->SetFilterMap( trk->GetFilterMap() );
 398       upcTrack->SetCharge( trk->Charge() );
 399       upcTrack->SetChi2perNDF( trk->Chi2perNDF() );
 400       upcTrack->SetTPCmomentum( apid->GetTPCmomentum() );
 401       upcTrack->SetTPCsignal( apid->GetTPCsignal() );
 402       upcTrack->SetTPCNcls( trk->GetTPCNcls() );
 403       upcTrack->SetTPCCrossedRows( (Float_t) trk->GetTPCNCrossedRows() );
 404       upcTrack->SetTPCNclsF( trk->GetTPCNclsF() );
 405       upcTrack->SetTPCNclsS( trk->GetTPCnclsS() );
 406       upcTrack->SetITSClusterMap( trk->GetITSClusterMap() );
 407
 408       //array for DCA
 409       //trackArray->Clear("C");
 410       //new((*trackArray)[0]) AliAODTrack(*trk);
 411       //new((*trackArray)[1]) AliAODTrack(*trk);
 412       //new((*trackArray)[2]) AliAODTrack(*trk);
 413
 414       //DCA to default primary vertex
 415       AliAODTrack *track1 = (AliAODTrack*) trk->Clone("track1");
 416       //AliAODTrack *track1 = (AliAODTrack*) trackArray->At( 0 );
 417       track1->PropagateToDCA(aodEvent->GetPrimaryVertex(), aodEvent->GetMagneticField(), 9999., xyzDca, cov);
 418       for(Int_t i=0; i<2; i++) {b[i] = (Float_t) xyzDca[i]; covF[i] = (Float_t) cov[i];}
 419       upcTrack->SetImpactParameters(b, covF);
 420       delete track1;
 421
 422       //DCA to SPD vertex
 423       track1 = (AliAODTrack*) trk->Clone("track1");
 424       //track1 = (AliAODTrack*) trackArray->At( 1 );
 425       track1->PropagateToDCA(aodEvent->GetPrimaryVertexSPD(), aodEvent->GetMagneticField(), 9999., xyzDca, cov);
 426       for(Int_t i=0; i<2; i++) {b[i] = (Float_t) xyzDca[i]; covF[i] = (Float_t) cov[i];}
 427       upcTrack->SetImpactParametersSPD(b, covF);
 428       delete track1;
 429
 430       //DCA to nominal interaction point
 431       track1 = (AliAODTrack*) trk->Clone("track1");
 432       //track1 = (AliAODTrack*) trackArray->At( 2 );
 433       track1->PropagateToDCA(vtx0, aodEvent->GetMagneticField(), 9999., xyzDca, cov);
 434       for(Int_t i=0; i<2; i++) {b[i] = (Float_t) xyzDca[i]; covF[i] = (Float_t) cov[i];}
 435       upcTrack->SetImpactParametersIP(b, covF);
 436       delete track1;
 437
 438       ncen++;
 439     }
 440   }// AOD tracks loop
 441
 442   //selection for at least one muon or central track
 443   if( nmun + ncen < 1 ) return kFALSE;
 444
 445   //selection for at least one muon and at least one central track
 446   //if( nmun < 1 || ncen < 1 ) return kFALSE;
 447
 448
 449   // L0 trigger inputs
 450   fUPCEvent->SetL0inputs( aodEvent->GetHeader()->GetL0TriggerInputs() );
 451
 452   // Tracklets
 453   fUPCEvent->SetNumberOfTracklets( aodEvent->GetTracklets()->GetNumberOfTracklets() );
 454
 455   // AOD ZDC for TDC
 456   AliAODZDC *dataZDCAOD = aodEvent->GetZDCData();
 457   if(!dataZDCAOD) {PostData(2, fHistList); return kFALSE;}
 458
 459   Bool_t znctdc = kFALSE, znatdc = kFALSE;
 460   if( dataZDCAOD->GetZNCTime() != 0. ) znctdc = kTRUE;
 461   if( dataZDCAOD->GetZNATime() != 0. ) znatdc = kTRUE;
 462   fUPCEvent->SetZNCtdc( znctdc );
 463   fUPCEvent->SetZNAtdc( znatdc );
 464   fUPCEvent->SetZPCtdc( kFALSE );
 465   fUPCEvent->SetZPAtdc( kFALSE );
 466
 467   //SPD primary vertex in AOD
 468   AliAODVertex *vtx = aodEvent->GetPrimaryVertexSPD();
 469   if(!vtx) {PostData(2, fHistList); return kFALSE;}
 470   Double_t posVtx[3]; vtx->GetXYZ( posVtx );
 471   Double_t covVtx[6]; vtx->GetCovarianceMatrix( covVtx );
 472   fUPCEvent->SetPrimaryVertexSPD(posVtx, vtx->GetChi2perNDF(), covVtx, vtx->GetNContributors());
 473   const char *vtxtitle = vtx->GetTitle();
 474   fUPCEvent->SetPrimaryVertexSPDtitle( vtxtitle );
 475
 476   return kTRUE;
 477
 478 }//RunAOD
 479
 480 //_____________________________________________________________________________
 481 void AliAnalysisTaskUpcFilterSemiforward::RunAODMC(TClonesArray *arrayMC, AliAODMCHeader *headerMC)
 482 {
 483   // run over AOD mc particles
 484
 485   if( !arrayMC || !headerMC ) return;
 486
 487   //generated vertex in AOD MC
 488   Double_t vtxD[3]; headerMC->GetVertex( vtxD );
 489   Float_t vtxF[3]; for(Int_t i=0; i<3; i++) vtxF[i] = (Float_t) vtxD[i];
 490   fUPCEvent->SetPrimaryVertexMC( vtxF );
 491
 492   //loop over mc particles
 493   for(Int_t imc=0; imc<arrayMC->GetEntriesFast(); imc++) {
 494     AliAODMCParticle *aodmc = (AliAODMCParticle*) arrayMC->At(imc);
 495     if(!aodmc) continue;
 496
 497     if(aodmc->GetMother() >= 0) continue;
 498
 499     TParticle *part = fUPCEvent->AddMCParticle();
 500     part->SetMomentum(aodmc->Px(), aodmc->Py(), aodmc->Pz(), aodmc->E());
 501     part->SetProductionVertex(aodmc->Xv(), aodmc->Yv(), aodmc->Zv(), aodmc->T());
 502     part->SetFirstMother(aodmc->GetMother());
 503     part->SetLastDaughter(aodmc->GetNDaughters());
 504     part->SetPdgCode(aodmc->GetPdgCode());
 505     part->SetUniqueID(imc);
 506   }//loop over mc particles
 507
 508 }//RunAODMC
 509
 510 //_____________________________________________________________________________
 511 Bool_t AliAnalysisTaskUpcFilterSemiforward::RunESD()
 512 {
 513   //cout<<"#################### ESD event ##################"<<endl;
 514
 515   // Input ESD event
 516   AliESDEvent *esdEvent = (AliESDEvent*) InputEvent();
 517   if(!esdEvent) return kFALSE;
 518
 519   fCounter->Fill( 21 ); // 21 = ESD analyzed events
 520
 521   if(fIsMC) {
 522     fUPCEvent->SetNSPDfiredInner( fTriggerAna->SPDFiredChips(esdEvent,1,kFALSE,1) );
 523     fUPCEvent->SetNSPDfiredOuter( fTriggerAna->SPDFiredChips(esdEvent,1,kFALSE,2) );
 524
 525     RunESDMC();
 526   }
 527
 528   static AliESDVertex *vtx0 = 0x0;
 529   if( !vtx0 ) {
 530     // make a vertex at coordinates 0, 0, 0
 531     vtx0 = new AliESDVertex();
 532     vtx0->SetXv(0.); vtx0->SetYv(0.); vtx0->SetZv(0.);
 533   }
 534
 535   // ESD central tracks
 536   Double_t pxpypz[3];
 537   Float_t b[2]; Float_t cov[3];
 538   UChar_t maskMan;
 539   UInt_t filterMap;
 540   Int_t nmun=0, ncen=0;
 541   //ESD central tracks loop
 542   for(Int_t itr=0; itr<esdEvent->GetNumberOfTracks(); itr++) {
 543     AliESDtrack *eTrack = esdEvent->GetTrack(itr);
 544     if( !eTrack ) continue;
 545
 546     // manual filter
 547     maskMan = 0;
 548     if( eTrack->HasPointOnITSLayer(0) )                 { maskMan |= 1 << 0; }
 549     if( eTrack->HasPointOnITSLayer(1) )                 { maskMan |= 1 << 1; }
 550     if( eTrack->GetStatus() & AliESDtrack::kITSrefit )  { maskMan |= 1 << 2; }
 551     if( eTrack->GetStatus() & AliESDtrack::kTPCrefit )  { maskMan |= 1 << 3; }
 552     if( eTrack->GetKinkIndex(0) > 0 )                   { maskMan |= 1 << 4; } // bit set if track is kink candidate
 553
 554     //selection for at least one point in ITS and TPC refit
 555     //if( !(maskMan & (1 << 0)) && !(maskMan & (1 << 1)) ) continue;
 556     //if( !(maskMan & (1 << 3)) ) continue;
 557
 558     //selection for its refit and tpc refit
 559     //if( !(maskMan & (1 << 2)) || !(maskMan & (1 << 3)) ) continue;
 560
 561     //selection for tpc refit
 562     //if( !(maskMan & (1 << 3)) ) continue;
 563
 564     //selection for at least one point in SPD and ITS refit and TPC refit
 565     //if( !(maskMan & (1 << 0)) && !(maskMan & (1 << 1)) ) continue;
 566     //if( !(maskMan & (1 << 2)) || !(maskMan & (1 << 3)) ) continue;
 567
 568     //selection for at least one point in SPD and at least one TPC cluster
 569     if( !(maskMan & (1 << 0)) && !(maskMan & (1 << 1)) ) continue;
 570     if( eTrack->GetTPCNcls() == 0 ) continue;
 571
 572     //central track accepted to write to the UPC event
 573
 574     // various configurations of AliESDtrackCuts
 575     filterMap = 0;
 576     for(UInt_t i=0; i<32; i++) {
 577       if( !fCutsList[i] ) continue;
 578
 579       if( fCutsList[i]->AcceptTrack(eTrack) ) { filterMap |= 1 << i; }
 580     }
 581
 582     eTrack->GetPxPyPz(pxpypz);
 583
 584     //fill UPC track
 585     AliUPCTrack *upcTrack = fUPCEvent->AddTrack();
 586     upcTrack->SetPxPyPz( pxpypz );
 587     upcTrack->SetMaskMan( maskMan );
 588     upcTrack->SetFilterMap( filterMap );
 589     upcTrack->SetCharge( eTrack->Charge() );
 590     upcTrack->SetChi2perNDF( eTrack->GetTPCchi2()/((Double_t) eTrack->GetTPCNcls()) ); // TPC chi2 used to fill this data member in esd
 591     upcTrack->SetTPCmomentum( eTrack->GetTPCmomentum() );
 592     upcTrack->SetTPCsignal( eTrack->GetTPCsignal() );
 593     upcTrack->SetTPCNcls( eTrack->GetTPCNcls() );
 594     upcTrack->SetTPCCrossedRows( eTrack->GetTPCCrossedRows() );
 595     upcTrack->SetTPCNclsF( eTrack->GetTPCNclsF() );
 596     upcTrack->SetTPCNclsS( eTrack->GetTPCnclsS() );
 597     upcTrack->SetITSchi2perNDF( eTrack->GetITSchi2()/((Double_t) eTrack->GetNcls(0)) );
 598     upcTrack->SetITSClusterMap( eTrack->GetITSClusterMap() );
 599
 600     //DCA to default primary vertex
 601     eTrack->GetImpactParameters(b, cov);
 602     upcTrack->SetImpactParameters(b, cov);
 603
 604     //DCA to SPD vertex
 605     AliESDtrack *track1 = (AliESDtrack*) eTrack->Clone("track1");
 606     track1->RelateToVertex( esdEvent->GetPrimaryVertexSPD(), esdEvent->GetMagneticField(), 9999. );
 607     track1->GetImpactParameters(b, cov);
 608     upcTrack->SetImpactParametersSPD(b, cov);
 609     delete track1;
 610
 611     //DCA to nominal interaction point
 612     track1 = (AliESDtrack*) eTrack->Clone("track1");
 613     track1->RelateToVertex( vtx0, esdEvent->GetMagneticField(), 9999. );
 614     track1->GetImpactParameters(b, cov);
 615     upcTrack->SetImpactParametersIP(b, cov);
 616     delete track1;
 617
 618     ncen++;
 619   } //ESD central tracks loop
 620
 621   // ESD muon tracks
 622   //muon tracks loop
 623   for(Int_t itr=0; itr<esdEvent->GetNumberOfMuonTracks(); itr++) {
 624     AliESDMuonTrack *mTrack = esdEvent->GetMuonTrack(itr);
 625     if( !mTrack ) continue;
 626
 627     AliUPCMuonTrack *upcMuon = fUPCEvent->AddMuonTrack();
 628     upcMuon->SetPtEtaPhi( mTrack->Pt(), mTrack->Eta(), mTrack->Phi() );
 629     upcMuon->SetCharge( mTrack->Charge() );
 630     upcMuon->SetMatchTrigger( mTrack->GetMatchTrigger() );
 631     upcMuon->SetRAtAbsorberEnd( mTrack->GetRAtAbsorberEnd() );
 632     upcMuon->SetChi2perNDF( mTrack->GetChi2()/((Double_t) mTrack->GetNDF()) );
 633     upcMuon->SetDCA( mTrack->GetDCA() );
 634     upcMuon->SetPxDCA( fMuonCuts->IsSelected(mTrack) );
 635
 636     nmun++;
 637   } //muon tracks loop
 638
 639   //selection for at least one muon or central track
 640   if( nmun + ncen < 1 ) return kFALSE;
 641
 642   //selection for at least one muon and at least one central track
 643   //if( nmun < 1 || ncen < 1 ) return kFALSE;
 644
 645   // L0 trigger inputs
 646   fUPCEvent->SetL0inputs( esdEvent->GetHeader()->GetL0TriggerInputs() );
 647
 648   //Tracklets
 649   fUPCEvent->SetNumberOfTracklets( esdEvent->GetMultiplicity()->GetNumberOfTracklets() );
 650
 651   // ESD ZDC for TDC
 652   AliESDZDC *dataZDCESD = esdEvent->GetESDZDC();
 653   if(!dataZDCESD) {PostData(2, fHistList); return kFALSE;}
 654
 655   Bool_t znctdc = kFALSE, zpctdc = kFALSE, znatdc = kFALSE, zpatdc = kFALSE;
 656   for(Int_t iz=0;iz<4;iz++) {
 657     if( dataZDCESD->GetZDCTDCData(10,iz) ) znctdc = kTRUE;
 658     if( dataZDCESD->GetZDCTDCData(11,iz) ) zpctdc = kTRUE;
 659     if( dataZDCESD->GetZDCTDCData(12,iz) ) znatdc = kTRUE;
 660     if( dataZDCESD->GetZDCTDCData(13,iz) ) zpatdc = kTRUE;
 661   }
 662   fUPCEvent->SetZNCtdc( znctdc );
 663   fUPCEvent->SetZPCtdc( zpctdc );
 664   fUPCEvent->SetZNAtdc( znatdc );
 665   fUPCEvent->SetZPAtdc( zpatdc );
 666
 667   //SPD primary vertex in ESD
 668   const AliESDVertex *vtx = esdEvent->GetPrimaryVertexSPD();
 669   if(!vtx) {PostData(2, fHistList); return kFALSE;}
 670   Double_t posVtx[3]; vtx->GetXYZ( posVtx );
 671   Double_t covVtx[6]; vtx->GetCovarianceMatrix( covVtx );
 672   fUPCEvent->SetPrimaryVertexSPD(posVtx, vtx->GetChi2perNDF(), covVtx, vtx->GetNContributors());
 673   const char *vtxtitle = vtx->GetTitle();
 674   fUPCEvent->SetPrimaryVertexSPDtitle( vtxtitle );
 675
 676   return kTRUE;
 677
 678 }//RunESD
 679
 680 //_____________________________________________________________________________
 681 void AliAnalysisTaskUpcFilterSemiforward::RunESDMC()
 682 {
 683   // ESD MC particles
 684
 685   AliMCEvent *mcEvent = MCEvent();
 686   if(!mcEvent) {cout<<"Error: no ESD MC found"<<endl; return;}
 687
 688   //generated vertex
 689   TArrayF vtx(3);
 690   mcEvent->GenEventHeader()->PrimaryVertex(vtx);
 691   fUPCEvent->SetPrimaryVertexMC( vtx.GetArray() );
 692
 693   //loop over mc particles
 694   for(Int_t imc=0; imc<mcEvent->GetNumberOfTracks(); imc++) {
 695     AliMCParticle *esdmc = (AliMCParticle*) mcEvent->GetTrack(imc);
 696     if(!esdmc) continue;
 697
 698     if(esdmc->GetMother() >= 0) continue;
 699
 700     TParticle *part = fUPCEvent->AddMCParticle();
 701     part->SetMomentum(esdmc->Px(), esdmc->Py(), esdmc->Pz(), esdmc->E());
 702     part->SetProductionVertex(esdmc->Xv(), esdmc->Yv(), esdmc->Zv(), 0.);
 703     part->SetFirstMother(esdmc->GetMother());
 704     part->SetLastDaughter(esdmc->GetLastDaughter()-esdmc->GetFirstDaughter()+1);
 705     part->SetPdgCode(esdmc->PdgCode());
 706     part->SetUniqueID(imc);
 707   }//loop over mc particles
 708
 709 }//RunESDMC
 710
 711 //_____________________________________________________________________________
 712 void AliAnalysisTaskUpcFilterSemiforward::Terminate(Option_t *)
 713 {
 714
 715   cout<<"Analysis complete."<<endl;
 716 }//Terminate
 717
 718
 719
 720
 721
 722
 723
 724
 725
 726
 727
 728
 729
 730
 731
 732
 733
 734
 735
 736
 737
 738
 739
 740
 741
 742
 743
 744
 745
 746