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", 44000, 154000, 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   trgClasses[11]= trigger.Contains("CMUP1-B"); // PbPb FW
 235
 236   Bool_t isTrg = kFALSE;
 237   for(Int_t itrg=1; itrg<NTRG; itrg++) {
 238     if(!trgClasses[itrg]) continue;
 239     //trigger at itrg is fired
 240     fUPCEvent->SetTriggerClass( itrg , kTRUE );
 241     fTriggerCounter->Fill( vEvent->GetRunNumber() , itrg );
 242     isTrg = kTRUE;
 243   }
 244   if(!isTrg && !fIsMC) {PostData(2, fHistList); return;}
 245   //event passed the trigger
 246
 247   fCounter->Fill( 2 ); // 2 = events after trigger (ESD and AOD)
 248
 249   // ESD / AOD specific tasks: MC, SPD FO, L0 inputs, SPD tracklets, tracks, ZDC tdc
 250   Bool_t stat;
 251   if( fIsESD ) stat = RunESD();
 252   if(!fIsESD ) stat = RunAOD();
 253   if( !stat && !fIsMC ) {PostData(2, fHistList); return;}
 254   //event passed ESD / AOD specific selection
 255
 256   fCounter->Fill( 3 ); // 3 = events after ESD / AOD specific part
 257
 258   // input data
 259   const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
 260   // reconstruction pass: -1 = unknown, 1 = pass1, 2 = pass2, counter indices: 9 (unknown), 11 (pass1), 12 (pass2)
 261   fUPCEvent->SetRecoPass( -1 );
 262   if( strstr(filnam,"/pass1/") ) {fUPCEvent->SetRecoPass( 1 ); fCounter->Fill( 11 );}
 263   if( strstr(filnam,"/pass2/") ) {fUPCEvent->SetRecoPass( 2 ); fCounter->Fill( 12 );}
 264   if( fUPCEvent->GetRecoPass() < 0 ) fCounter->Fill( 9 );
 265   fUPCEvent->SetInputFileName( filnam );
 266   fUPCEvent->SetEventNumber( ((TTree*) GetInputData(0))->GetTree()->GetReadEntry() );
 267   fUPCEvent->SetRunNumber( vEvent->GetRunNumber() );
 268
 269   //VZERO
 270   AliVVZERO *dataVZERO = vEvent->GetVZEROData();
 271   if(!dataVZERO) {PostData(2, fHistList); return;}
 272
 273   fUPCEvent->SetV0ADecision( dataVZERO->GetV0ADecision() );
 274   fUPCEvent->SetV0CDecision( dataVZERO->GetV0CDecision() );
 275   for (UInt_t iv=0; iv<32; iv++) {
 276     if( dataVZERO->BBTriggerV0C((Int_t)iv) ) fUPCEvent->SetBBtriggerV0Cmask(iv);
 277     if( dataVZERO->GetBBFlag((Int_t)iv) ) fUPCEvent->SetBBFlagV0Cmask(iv);
 278   }
 279
 280   //ZDC
 281   AliVZDC *dataZDC = vEvent->GetZDCData();
 282   if(!dataZDC) {PostData(2, fHistList); return;}
 283
 284   //energy in ZDC
 285   Double_t eZnc=0., eZpc=0., eZna=0., eZpa=0.;
 286   for(Int_t i=0; i<5; i++) {
 287     eZnc += dataZDC->GetZNCTowerEnergy()[i];
 288     eZpc += dataZDC->GetZPCTowerEnergy()[i];
 289     eZna += dataZDC->GetZNATowerEnergy()[i];
 290     eZpa += dataZDC->GetZPATowerEnergy()[i];
 291   }
 292   fUPCEvent->SetZNCEnergy( eZnc );
 293   fUPCEvent->SetZPCEnergy( eZpc );
 294   fUPCEvent->SetZNAEnergy( eZna );
 295   fUPCEvent->SetZPAEnergy( eZpa );
 296
 297   //default primary vertex
 298   const AliVVertex *vtx = vEvent->GetPrimaryVertex();
 299   if(!vtx) {PostData(2, fHistList); return;}
 300   Double_t pos[3]; vtx->GetXYZ(pos);
 301   Double_t cov[6]; vtx->GetCovarianceMatrix(cov);
 302   fUPCEvent->SetPrimaryVertex(pos, vtx->GetChi2perNDF(), cov, vtx->GetNContributors());
 303   const char *vtxtitle = vtx->GetTitle();
 304   fUPCEvent->SetPrimaryVertexTitle( vtxtitle );
 305
 306   fCounter->Fill( 4 ); // 4 = events written to the tree (ESD and AOD)
 307
 308   fUPCTree ->Fill();
 309   PostData(1, fUPCTree);
 310   PostData(2, fHistList);
 311
 312 }//UserExec
 313
 314 //_____________________________________________________________________________
 315 Bool_t AliAnalysisTaskUpcFilterSemiforward::RunAOD()
 316 {
 317   //cout<<"#################### AOD event ##################"<<endl;
 318
 319   //input AOD event
 320   AliAODEvent *aodEvent = (AliAODEvent*) InputEvent();
 321   if(!aodEvent) return kFALSE;
 322
 323   fCounter->Fill( 22 ); // 22 = AOD analyzed events
 324
 325   if(fIsMC) RunAODMC( (TClonesArray*) aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName()),
 326                       (AliAODMCHeader*) aodEvent->FindListObject(AliAODMCHeader::StdBranchName())
 327                     );
 328
 329   //nominal interaction point
 330   static AliAODVertex *vtx0 = 0x0;
 331   if( !vtx0 ) {
 332     // make a vertex at coordinates 0, 0, 0
 333     vtx0 = new AliAODVertex();
 334     vtx0->SetX(0.); vtx0->SetY(0.); vtx0->SetZ(0.);
 335   }
 336
 337   //array of tracks for DCAs calculation
 338   //static TClonesArray *trackArray = 0x0;
 339   //if( !trackArray ) trackArray = new TClonesArray("AliAODTrack", 3); // number of DCAs
 340
 341   Double_t pxpypz[3];
 342   UChar_t maskMan;
 343   Double_t xyzDca[2], cov[3];
 344   Float_t b[2], covF[3];
 345   Int_t nmun=0, ncen=0;
 346   // AOD tracks loop
 347   for(Int_t itr=0; itr<aodEvent->GetNumberOfTracks(); itr++) {
 348     AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aodEvent->GetTrack(itr));
 349     if( !trk ) continue;
 350
 351     //muon track
 352     if( trk->IsMuonTrack() ) {
 353
 354       AliUPCMuonTrack *upcMuon = fUPCEvent->AddMuonTrack();
 355       upcMuon->SetPtEtaPhi( trk->Pt(), trk->Eta(), trk->Phi() );
 356       upcMuon->SetCharge( trk->Charge() );
 357       upcMuon->SetMatchTrigger( trk->GetMatchTrigger() );
 358       upcMuon->SetRAtAbsorberEnd( trk->GetRAtAbsorberEnd() );
 359       upcMuon->SetChi2perNDF( trk->Chi2perNDF() );
 360       upcMuon->SetDCA( trk->DCA() );
 361       upcMuon->SetPxDCA( fMuonCuts->IsSelected(trk) );
 362
 363       nmun++;
 364     }
 365     //central track
 366     else {
 367
 368       maskMan = 0;
 369       if( trk->HasPointOnITSLayer(0) )                 { maskMan |= 1 << 0; }
 370       if( trk->HasPointOnITSLayer(1) )                 { maskMan |= 1 << 1; }
 371       if( trk->GetStatus() & AliESDtrack::kITSrefit )  { maskMan |= 1 << 2; }
 372       if( trk->GetStatus() & AliESDtrack::kTPCrefit )  { maskMan |= 1 << 3; }
 373
 374       //selection for at least one point in ITS and TPC refit
 375       //if( !(maskMan & (1 << 0)) && !(maskMan & (1 << 1)) ) continue;
 376       //if( !(maskMan & (1 << 3)) ) continue;
 377
 378       //selection for its refit and tpc refit
 379       //if( !(maskMan & (1 << 2)) || !(maskMan & (1 << 3)) ) continue;
 380
 381       //selection for tpc refit
 382       //if( !(maskMan & (1 << 3)) ) continue;
 383
 384       //selection for at least one point in SPD and ITS refit and TPC refit
 385       //if( !(maskMan & (1 << 0)) && !(maskMan & (1 << 1)) ) continue;
 386       //if( !(maskMan & (1 << 2)) || !(maskMan & (1 << 3)) ) continue;
 387
 388       //selection for at least one point in SPD and at least one TPC cluster
 389       if( !(maskMan & (1 << 0)) && !(maskMan & (1 << 1)) ) continue;
 390       if( trk->GetTPCNcls() == 0 ) continue;
 391
 392       trk->GetPxPyPz(pxpypz);
 393       AliAODPid *apid = trk->GetDetPid();
 394       if(!apid) continue;
 395
 396       AliUPCTrack *upcTrack = fUPCEvent->AddTrack();
 397       upcTrack->SetPxPyPz( pxpypz );
 398       upcTrack->SetMaskMan( maskMan );
 399       upcTrack->SetFilterMap( trk->GetFilterMap() );
 400       upcTrack->SetCharge( trk->Charge() );
 401       upcTrack->SetChi2perNDF( trk->Chi2perNDF() );
 402       upcTrack->SetTPCmomentum( apid->GetTPCmomentum() );
 403       upcTrack->SetTPCsignal( apid->GetTPCsignal() );
 404       upcTrack->SetTPCNcls( trk->GetTPCNcls() );
 405       upcTrack->SetTPCCrossedRows( (Float_t) trk->GetTPCNCrossedRows() );
 406       upcTrack->SetTPCNclsF( trk->GetTPCNclsF() );
 407       upcTrack->SetTPCNclsS( trk->GetTPCnclsS() );
 408       upcTrack->SetITSClusterMap( trk->GetITSClusterMap() );
 409
 410       //array for DCA
 411       //trackArray->Clear("C");
 412       //new((*trackArray)[0]) AliAODTrack(*trk);
 413       //new((*trackArray)[1]) AliAODTrack(*trk);
 414       //new((*trackArray)[2]) AliAODTrack(*trk);
 415
 416       //DCA to default primary vertex
 417       AliAODTrack *track1 = (AliAODTrack*) trk->Clone("track1");
 418       //AliAODTrack *track1 = (AliAODTrack*) trackArray->At( 0 );
 419       track1->PropagateToDCA(aodEvent->GetPrimaryVertex(), aodEvent->GetMagneticField(), 9999., xyzDca, cov);
 420       for(Int_t i=0; i<2; i++) {b[i] = (Float_t) xyzDca[i]; covF[i] = (Float_t) cov[i];}
 421       upcTrack->SetImpactParameters(b, covF);
 422       delete track1;
 423
 424       //DCA to SPD vertex
 425       track1 = (AliAODTrack*) trk->Clone("track1");
 426       //track1 = (AliAODTrack*) trackArray->At( 1 );
 427       track1->PropagateToDCA(aodEvent->GetPrimaryVertexSPD(), aodEvent->GetMagneticField(), 9999., xyzDca, cov);
 428       for(Int_t i=0; i<2; i++) {b[i] = (Float_t) xyzDca[i]; covF[i] = (Float_t) cov[i];}
 429       upcTrack->SetImpactParametersSPD(b, covF);
 430       delete track1;
 431
 432       //DCA to nominal interaction point
 433       track1 = (AliAODTrack*) trk->Clone("track1");
 434       //track1 = (AliAODTrack*) trackArray->At( 2 );
 435       track1->PropagateToDCA(vtx0, aodEvent->GetMagneticField(), 9999., xyzDca, cov);
 436       for(Int_t i=0; i<2; i++) {b[i] = (Float_t) xyzDca[i]; covF[i] = (Float_t) cov[i];}
 437       upcTrack->SetImpactParametersIP(b, covF);
 438       delete track1;
 439
 440       ncen++;
 441     }
 442   }// AOD tracks loop
 443
 444   //selection for at least one muon or central track
 445   if( nmun + ncen < 1 ) return kFALSE;
 446
 447   //selection for at least one muon and at least one central track
 448   //if( nmun < 1 || ncen < 1 ) return kFALSE;
 449
 450
 451   // L0 trigger inputs
 452   fUPCEvent->SetL0inputs( aodEvent->GetHeader()->GetL0TriggerInputs() );
 453
 454   // Tracklets
 455   fUPCEvent->SetNumberOfTracklets( aodEvent->GetTracklets()->GetNumberOfTracklets() );
 456
 457   // AOD ZDC for TDC
 458   AliAODZDC *dataZDCAOD = aodEvent->GetZDCData();
 459   if(!dataZDCAOD) {PostData(2, fHistList); return kFALSE;}
 460
 461   Bool_t znctdc = kFALSE, znatdc = kFALSE;
 462   if( dataZDCAOD->GetZNCTime() != 0. ) znctdc = kTRUE;
 463   if( dataZDCAOD->GetZNATime() != 0. ) znatdc = kTRUE;
 464   fUPCEvent->SetZNCtdc( znctdc );
 465   fUPCEvent->SetZNAtdc( znatdc );
 466   fUPCEvent->SetZPCtdc( kFALSE );
 467   fUPCEvent->SetZPAtdc( kFALSE );
 468
 469   //SPD primary vertex in AOD
 470   AliAODVertex *vtx = aodEvent->GetPrimaryVertexSPD();
 471   if(!vtx) {PostData(2, fHistList); return kFALSE;}
 472   Double_t posVtx[3]; vtx->GetXYZ( posVtx );
 473   Double_t covVtx[6]; vtx->GetCovarianceMatrix( covVtx );
 474   fUPCEvent->SetPrimaryVertexSPD(posVtx, vtx->GetChi2perNDF(), covVtx, vtx->GetNContributors());
 475   const char *vtxtitle = vtx->GetTitle();
 476   fUPCEvent->SetPrimaryVertexSPDtitle( vtxtitle );
 477
 478   return kTRUE;
 479
 480 }//RunAOD
 481
 482 //_____________________________________________________________________________
 483 void AliAnalysisTaskUpcFilterSemiforward::RunAODMC(TClonesArray *arrayMC, AliAODMCHeader *headerMC)
 484 {
 485   // run over AOD mc particles
 486
 487   if( !arrayMC || !headerMC ) return;
 488
 489   //generated vertex in AOD MC
 490   Double_t vtxD[3]; headerMC->GetVertex( vtxD );
 491   Float_t vtxF[3]; for(Int_t i=0; i<3; i++) vtxF[i] = (Float_t) vtxD[i];
 492   fUPCEvent->SetPrimaryVertexMC( vtxF );
 493
 494   //loop over mc particles
 495   for(Int_t imc=0; imc<arrayMC->GetEntriesFast(); imc++) {
 496     AliAODMCParticle *aodmc = (AliAODMCParticle*) arrayMC->At(imc);
 497     if(!aodmc) continue;
 498
 499     if(aodmc->GetMother() >= 0) continue;
 500
 501     TParticle *part = fUPCEvent->AddMCParticle();
 502     part->SetMomentum(aodmc->Px(), aodmc->Py(), aodmc->Pz(), aodmc->E());
 503     part->SetProductionVertex(aodmc->Xv(), aodmc->Yv(), aodmc->Zv(), aodmc->T());
 504     part->SetFirstMother(aodmc->GetMother());
 505     part->SetLastDaughter(aodmc->GetNDaughters());
 506     part->SetPdgCode(aodmc->GetPdgCode());
 507     part->SetUniqueID(imc);
 508   }//loop over mc particles
 509
 510 }//RunAODMC
 511
 512 //_____________________________________________________________________________
 513 Bool_t AliAnalysisTaskUpcFilterSemiforward::RunESD()
 514 {
 515   //cout<<"#################### ESD event ##################"<<endl;
 516
 517   // Input ESD event
 518   AliESDEvent *esdEvent = (AliESDEvent*) InputEvent();
 519   if(!esdEvent) return kFALSE;
 520
 521   fCounter->Fill( 21 ); // 21 = ESD analyzed events
 522
 523   if(fIsMC) {
 524     fUPCEvent->SetNSPDfiredInner( fTriggerAna->SPDFiredChips(esdEvent,1,kFALSE,1) );
 525     fUPCEvent->SetNSPDfiredOuter( fTriggerAna->SPDFiredChips(esdEvent,1,kFALSE,2) );
 526
 527     RunESDMC();
 528   }
 529
 530   static AliESDVertex *vtx0 = 0x0;
 531   if( !vtx0 ) {
 532     // make a vertex at coordinates 0, 0, 0
 533     vtx0 = new AliESDVertex();
 534     vtx0->SetXv(0.); vtx0->SetYv(0.); vtx0->SetZv(0.);
 535   }
 536
 537   // ESD central tracks
 538   Double_t pxpypz[3];
 539   Float_t b[2]; Float_t cov[3];
 540   UChar_t maskMan;
 541   UInt_t filterMap;
 542   Int_t nmun=0, ncen=0;
 543   //ESD central tracks loop
 544   for(Int_t itr=0; itr<esdEvent->GetNumberOfTracks(); itr++) {
 545     AliESDtrack *eTrack = esdEvent->GetTrack(itr);
 546     if( !eTrack ) continue;
 547
 548     // manual filter
 549     maskMan = 0;
 550     if( eTrack->HasPointOnITSLayer(0) )                 { maskMan |= 1 << 0; }
 551     if( eTrack->HasPointOnITSLayer(1) )                 { maskMan |= 1 << 1; }
 552     if( eTrack->GetStatus() & AliESDtrack::kITSrefit )  { maskMan |= 1 << 2; }
 553     if( eTrack->GetStatus() & AliESDtrack::kTPCrefit )  { maskMan |= 1 << 3; }
 554     if( eTrack->GetKinkIndex(0) > 0 )                   { maskMan |= 1 << 4; } // bit set if track is kink candidate
 555
 556     //selection for at least one point in ITS and TPC refit
 557     //if( !(maskMan & (1 << 0)) && !(maskMan & (1 << 1)) ) continue;
 558     //if( !(maskMan & (1 << 3)) ) continue;
 559
 560     //selection for its refit and tpc refit
 561     //if( !(maskMan & (1 << 2)) || !(maskMan & (1 << 3)) ) continue;
 562
 563     //selection for tpc refit
 564     //if( !(maskMan & (1 << 3)) ) continue;
 565
 566     //selection for at least one point in SPD and ITS refit and TPC refit
 567     //if( !(maskMan & (1 << 0)) && !(maskMan & (1 << 1)) ) continue;
 568     //if( !(maskMan & (1 << 2)) || !(maskMan & (1 << 3)) ) continue;
 569
 570     //selection for at least one point in SPD and at least one TPC cluster
 571     if( !(maskMan & (1 << 0)) && !(maskMan & (1 << 1)) ) continue;
 572     if( eTrack->GetTPCNcls() == 0 ) continue;
 573
 574     //central track accepted to write to the UPC event
 575
 576     // various configurations of AliESDtrackCuts
 577     filterMap = 0;
 578     for(UInt_t i=0; i<32; i++) {
 579       if( !fCutsList[i] ) continue;
 580
 581       if( fCutsList[i]->AcceptTrack(eTrack) ) { filterMap |= 1 << i; }
 582     }
 583
 584     eTrack->GetPxPyPz(pxpypz);
 585
 586     //fill UPC track
 587     AliUPCTrack *upcTrack = fUPCEvent->AddTrack();
 588     upcTrack->SetPxPyPz( pxpypz );
 589     upcTrack->SetMaskMan( maskMan );
 590     upcTrack->SetFilterMap( filterMap );
 591     upcTrack->SetCharge( eTrack->Charge() );
 592     upcTrack->SetChi2perNDF( eTrack->GetTPCchi2()/((Double_t) eTrack->GetTPCNcls()) ); // TPC chi2 used to fill this data member in esd
 593     upcTrack->SetTPCmomentum( eTrack->GetTPCmomentum() );
 594     upcTrack->SetTPCsignal( eTrack->GetTPCsignal() );
 595     upcTrack->SetTPCNcls( eTrack->GetTPCNcls() );
 596     upcTrack->SetTPCCrossedRows( eTrack->GetTPCCrossedRows() );
 597     upcTrack->SetTPCNclsF( eTrack->GetTPCNclsF() );
 598     upcTrack->SetTPCNclsS( eTrack->GetTPCnclsS() );
 599     upcTrack->SetITSchi2perNDF( eTrack->GetITSchi2()/((Double_t) eTrack->GetNcls(0)) );
 600     upcTrack->SetITSClusterMap( eTrack->GetITSClusterMap() );
 601
 602     //DCA to default primary vertex
 603     eTrack->GetImpactParameters(b, cov);
 604     upcTrack->SetImpactParameters(b, cov);
 605
 606     //DCA to SPD vertex
 607     AliESDtrack *track1 = (AliESDtrack*) eTrack->Clone("track1");
 608     track1->RelateToVertex( esdEvent->GetPrimaryVertexSPD(), esdEvent->GetMagneticField(), 9999. );
 609     track1->GetImpactParameters(b, cov);
 610     upcTrack->SetImpactParametersSPD(b, cov);
 611     delete track1;
 612
 613     //DCA to nominal interaction point
 614     track1 = (AliESDtrack*) eTrack->Clone("track1");
 615     track1->RelateToVertex( vtx0, esdEvent->GetMagneticField(), 9999. );
 616     track1->GetImpactParameters(b, cov);
 617     upcTrack->SetImpactParametersIP(b, cov);
 618     delete track1;
 619
 620     ncen++;
 621   } //ESD central tracks loop
 622
 623   // ESD muon tracks
 624   //muon tracks loop
 625   for(Int_t itr=0; itr<esdEvent->GetNumberOfMuonTracks(); itr++) {
 626     AliESDMuonTrack *mTrack = esdEvent->GetMuonTrack(itr);
 627     if( !mTrack ) continue;
 628
 629     AliUPCMuonTrack *upcMuon = fUPCEvent->AddMuonTrack();
 630     upcMuon->SetPtEtaPhi( mTrack->Pt(), mTrack->Eta(), mTrack->Phi() );
 631     upcMuon->SetCharge( mTrack->Charge() );
 632     upcMuon->SetMatchTrigger( mTrack->GetMatchTrigger() );
 633     upcMuon->SetRAtAbsorberEnd( mTrack->GetRAtAbsorberEnd() );
 634     upcMuon->SetChi2perNDF( mTrack->GetChi2()/((Double_t) mTrack->GetNDF()) );
 635     upcMuon->SetDCA( mTrack->GetDCA() );
 636     upcMuon->SetPxDCA( fMuonCuts->IsSelected(mTrack) );
 637
 638     nmun++;
 639   } //muon tracks loop
 640
 641   //selection for at least one muon or central track
 642   if( nmun + ncen < 1 ) return kFALSE;
 643
 644   //selection for at least one muon and at least one central track
 645   //if( nmun < 1 || ncen < 1 ) return kFALSE;
 646
 647   // L0 trigger inputs
 648   fUPCEvent->SetL0inputs( esdEvent->GetHeader()->GetL0TriggerInputs() );
 649
 650   //Tracklets
 651   fUPCEvent->SetNumberOfTracklets( esdEvent->GetMultiplicity()->GetNumberOfTracklets() );
 652
 653   // ESD ZDC for TDC
 654   AliESDZDC *dataZDCESD = esdEvent->GetESDZDC();
 655   if(!dataZDCESD) {PostData(2, fHistList); return kFALSE;}
 656
 657   Bool_t znctdc = kFALSE, zpctdc = kFALSE, znatdc = kFALSE, zpatdc = kFALSE;
 658   for(Int_t iz=0;iz<4;iz++) {
 659     if( dataZDCESD->GetZDCTDCData(10,iz) ) znctdc = kTRUE;
 660     if( dataZDCESD->GetZDCTDCData(11,iz) ) zpctdc = kTRUE;
 661     if( dataZDCESD->GetZDCTDCData(12,iz) ) znatdc = kTRUE;
 662     if( dataZDCESD->GetZDCTDCData(13,iz) ) zpatdc = kTRUE;
 663   }
 664   fUPCEvent->SetZNCtdc( znctdc );
 665   fUPCEvent->SetZPCtdc( zpctdc );
 666   fUPCEvent->SetZNAtdc( znatdc );
 667   fUPCEvent->SetZPAtdc( zpatdc );
 668
 669   //SPD primary vertex in ESD
 670   const AliESDVertex *vtx = esdEvent->GetPrimaryVertexSPD();
 671   if(!vtx) {PostData(2, fHistList); return kFALSE;}
 672   Double_t posVtx[3]; vtx->GetXYZ( posVtx );
 673   Double_t covVtx[6]; vtx->GetCovarianceMatrix( covVtx );
 674   fUPCEvent->SetPrimaryVertexSPD(posVtx, vtx->GetChi2perNDF(), covVtx, vtx->GetNContributors());
 675   const char *vtxtitle = vtx->GetTitle();
 676   fUPCEvent->SetPrimaryVertexSPDtitle( vtxtitle );
 677
 678   return kTRUE;
 679
 680 }//RunESD
 681
 682 //_____________________________________________________________________________
 683 void AliAnalysisTaskUpcFilterSemiforward::RunESDMC()
 684 {
 685   // ESD MC particles
 686
 687   AliMCEvent *mcEvent = MCEvent();
 688   if(!mcEvent) {cout<<"Error: no ESD MC found"<<endl; return;}
 689
 690   //generated vertex
 691   TArrayF vtx(3);
 692   mcEvent->GenEventHeader()->PrimaryVertex(vtx);
 693   fUPCEvent->SetPrimaryVertexMC( vtx.GetArray() );
 694
 695   //loop over mc particles
 696   for(Int_t imc=0; imc<mcEvent->GetNumberOfTracks(); imc++) {
 697     AliMCParticle *esdmc = (AliMCParticle*) mcEvent->GetTrack(imc);
 698     if(!esdmc) continue;
 699
 700     if(esdmc->GetMother() >= 0) continue;
 701
 702     TParticle *part = fUPCEvent->AddMCParticle();
 703     part->SetMomentum(esdmc->Px(), esdmc->Py(), esdmc->Pz(), esdmc->E());
 704     part->SetProductionVertex(esdmc->Xv(), esdmc->Yv(), esdmc->Zv(), 0.);
 705     part->SetFirstMother(esdmc->GetMother());
 706     part->SetLastDaughter(esdmc->GetLastDaughter()-esdmc->GetFirstDaughter()+1);
 707     part->SetPdgCode(esdmc->PdgCode());
 708     part->SetUniqueID(imc);
 709   }//loop over mc particles
 710
 711 }//RunESDMC
 712
 713 //_____________________________________________________________________________
 714 void AliAnalysisTaskUpcFilterSemiforward::Terminate(Option_t *)
 715 {
 716
 717   cout<<"Analysis complete."<<endl;
 718 }//Terminate
 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
 747
 748