PWGUD/UPC/AliAnalysisTaskUpcK0sK0s.cxx

   1 /*************************************************************************
   2 * Copyright(c) 1998-2008, 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 // c++ headers
  17 #include <iostream>
  18 #include <string.h>
  19 #include <bitset>
  20
  21 // root headers
  22 #include "TH1I.h"
  23 #include "TTree.h"
  24 #include "TClonesArray.h"
  25 #include "TParticle.h"
  26 #include "TString.h"
  27 #include "TFile.h"
  28 #include "TDatabasePDG.h"
  29 #include "TLorentzVector.h"
  30
  31 // aliroot headers
  32 #include "AliAnalysisManager.h"
  33 #include "AliInputEventHandler.h"
  34 #include "AliESDEvent.h"
  35 #include "AliAODEvent.h"
  36 #include "AliMCEvent.h"
  37 #include "AliAODVZERO.h"
  38 #include "AliAODZDC.h"
  39 #include "AliESDVZERO.h"
  40 #include "AliESDZDC.h"
  41 #include "AliPIDResponse.h"
  42 #include "AliAODTrack.h"
  43 #include "AliAODPid.h"
  44 #include "AliAODVertex.h"
  45 #include "AliESDVertex.h"
  46 #include "AliMultiplicity.h"
  47 #include "AliESDtrack.h"
  48 #include "AliESDMuonTrack.h"
  49 #include "AliAODMCParticle.h"
  50 #include "AliMCParticle.h"
  51
  52 // my headers
  53 #include "AliAnalysisTaskUpcK0sK0s.h"
  54
  55 ClassImp(AliAnalysisTaskUpcK0sK0s);
  56
  57 using std::cout;
  58 using std::endl;
  59
  60 //trees for UPC analysis,
  61 // michal.broz@cern.ch
  62
  63 //_____________________________________________________________________________
  64 AliAnalysisTaskUpcK0sK0s::AliAnalysisTaskUpcK0sK0s()
  65   : AliAnalysisTaskSE(),fType(0),fRunTree(kTRUE),fRunHist(kTRUE),fK0sTree(0),
  66     fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),fVtxContrib(0),fBCrossNum(0),fNtracklets(0),
  67     fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
  68     fDataFilnam(0),fRecoPass(0),fEvtNum(0),
  69     fK0sAODv0s(0),fK0sESDv0s(0),fK0sAODTracks(0),fK0sESDTracks(0),
  70     fListHist(0),fHistTriggersPerRun(0),fHistK0sMass(0)
  71
  72 {
  73
  74 //Dummy constructor
  75
  76 }//AliAnalysisTaskUpcK0sK0s
  77
  78
  79 //_____________________________________________________________________________
  80 AliAnalysisTaskUpcK0sK0s::AliAnalysisTaskUpcK0sK0s(const char *name)
  81   : AliAnalysisTaskSE(name),fType(0),fRunTree(kTRUE),fRunHist(kTRUE),fK0sTree(0),
  82     fRunNum(0),fPerNum(0),fOrbNum(0),fL0inputs(0),fL1inputs(0),fVtxContrib(0),fBCrossNum(0),fNtracklets(0),
  83     fZDCAenergy(0),fZDCCenergy(0),fV0Adecision(0),fV0Cdecision(0),
  84     fDataFilnam(0),fRecoPass(0),fEvtNum(0),
  85     fK0sAODv0s(0),fK0sESDv0s(0),fK0sAODTracks(0),fK0sESDTracks(0),
  86     fListHist(0),fHistTriggersPerRun(0),fHistK0sMass(0)
  87
  88 {
  89
  90   // Constructor
  91   if( strstr(name,"ESD") ) fType = 0;
  92   if( strstr(name,"AOD") ) fType = 1;
  93
  94   Init();
  95
  96   DefineOutput(1, TTree::Class());
  97   DefineOutput(2, TList::Class());
  98
  99 }//AliAnalysisTaskUpcK0sK0s
 100
 101 //_____________________________________________________________________________
 102 void AliAnalysisTaskUpcK0sK0s::Init()
 103 {
 104
 105   for(Int_t i=0; i<ntrg; i++) fTrigger[i] = kFALSE;
 106
 107 }//Init
 108
 109 //_____________________________________________________________________________
 110 AliAnalysisTaskUpcK0sK0s::~AliAnalysisTaskUpcK0sK0s()
 111 {
 112   // Destructor
 113   if(fK0sTree){
 114      delete fK0sTree;
 115      fK0sTree = 0x0;
 116   }
 117   if(fListHist){
 118      delete fListHist;
 119      fListHist = 0x0;
 120   }
 121
 122 }//~AliAnalysisTaskUpcK0sK0s
 123
 124
 125 //_____________________________________________________________________________
 126 void AliAnalysisTaskUpcK0sK0s::UserCreateOutputObjects()
 127 {
 128
 129     //vertices
 130   fK0sAODv0s = new TClonesArray("AliAODv0", 1000);
 131   fK0sESDv0s = new TClonesArray("AliESDv0", 1000);
 132
 133     //tracks
 134   fK0sAODTracks = new TClonesArray("AliAODTrack", 1000);
 135   fK0sESDTracks = new TClonesArray("AliESDtrack", 1000);
 136
 137   //output tree with K0s candidate events
 138   fK0sTree = new TTree("fK0sTree", "fK0sTree");
 139   fK0sTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
 140   fK0sTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
 141   fK0sTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
 142
 143   fK0sTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
 144   fK0sTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
 145   fK0sTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
 146   fK0sTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
 147   fK0sTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
 148   fK0sTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
 149   fK0sTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
 150   fK0sTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
 151   fK0sTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
 152   fK0sTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");
 153   fK0sTree ->Branch("fDataFilnam", &fDataFilnam);
 154   fK0sTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
 155   fK0sTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");
 156   if( fType == 0 ) {
 157     fK0sTree ->Branch("fK0sESDv0s", &fK0sESDv0s);
 158     fK0sTree ->Branch("fK0sESDTracks", &fK0sESDTracks);
 159   }
 160   if( fType == 1 ) {
 161     fK0sTree ->Branch("fK0sAODv0s", &fK0sAODv0s);
 162     fK0sTree ->Branch("fK0sAODTracks", &fK0sAODTracks);
 163   }
 164
 165
 166   fListHist = new TList();
 167   fListHist ->SetOwner();
 168
 169   fHistTriggersPerRun = new TH1D("fHistTriggersPerRun", "fHistTriggersPerRun", 3000, 167000.5, 170000.5);
 170   fListHist->Add(fHistTriggersPerRun);
 171
 172   fHistK0sMass = new TH2D("fHistK0sMass","fHistK0sMass",20,0.45,0.55,20,0.45,0.55);
 173   fListHist->Add(fHistK0sMass);
 174
 175   PostData(1, fK0sTree);
 176   PostData(2, fListHist);
 177
 178 }//UserCreateOutputObjects
 179
 180 //_____________________________________________________________________________
 181 void AliAnalysisTaskUpcK0sK0s::UserExec(Option_t *)
 182 {
 183
 184   //cout<<"#################### Next event ##################"<<endl;
 185
 186   if( fType == 0 ) RunESD();
 187   if( fType == 1 ){
 188         if(fRunHist) RunAODhist();
 189         if(fRunTree) RunAODtree();
 190         }
 191
 192 }//UserExec
 193 //_____________________________________________________________________________
 194 void AliAnalysisTaskUpcK0sK0s::RunAODhist()
 195 {
 196
 197   //input event
 198   AliAODEvent *aod = (AliAODEvent*) InputEvent();
 199   if(!aod) return;
 200
 201   //Trigger
 202   TString trigger = aod->GetFiredTriggerClasses();
 203
 204   if( !trigger.Contains("CCUP") ) return;
 205
 206   fRunNum = aod ->GetRunNumber();
 207   fHistTriggersPerRun->Fill(fRunNum);
 208
 209   //primary vertex
 210   AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
 211   fVtxContrib = fAODVertex->GetNContributors();
 212   if(fVtxContrib < 2) return;
 213
 214
 215   //VZERO, ZDC
 216   AliAODVZERO *fV0data = aod ->GetVZEROData();
 217
 218   fV0Adecision = fV0data->GetV0ADecision();
 219   fV0Cdecision = fV0data->GetV0CDecision();
 220   if(fV0Adecision != AliAODVZERO::kV0Empty || fV0Cdecision != AliAODVZERO::kV0Empty) return;
 221
 222   AliAODZDC *fZDCdata = aod->GetZDCData();
 223
 224   fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
 225   fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
 226   if( fZDCAenergy > 8200 || fZDCCenergy > 8200) return;
 227
 228   Int_t nGoodV0s=0;
 229   Int_t V0Index[3] = {-1,-1,-1};
 230   Int_t V0TrackID[6] = {-1,-1,-1,-1,-1,-1};
 231
 232   Int_t nGoodTracks=0;
 233   Int_t TrackID[5] = {-1,-1,-1,-1,-1};
 234
 235
 236   //V0s loop
 237   for(Int_t iV0=0; iV0<aod ->GetNumberOfV0s(); iV0++) {
 238     AliAODv0 *v0 = aod->GetV0(iV0);
 239     if( !v0 ) continue;
 240     Bool_t lOnFlyStatus = v0->GetOnFlyStatus();
 241     if (lOnFlyStatus) continue;
 242
 243     AliAODTrack *pTrack=(AliAODTrack *)v0->GetDaughter(0); //0->Positive Daughter
 244     AliAODTrack *nTrack=(AliAODTrack *)v0->GetDaughter(1); //1->Negative Daughter
 245     if (!pTrack || !nTrack) continue;
 246
 247     if ( pTrack->Charge() == nTrack->Charge())continue;
 248
 249       if(!(pTrack->TestFilterBit(1<<0))) continue;
 250       if(!(nTrack->TestFilterBit(1<<0))) continue;
 251       if(!(pTrack->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 252       if(!(nTrack->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 253       if(!(pTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 254       if(!(nTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 255       if(pTrack->GetTPCNcls() < 50)continue;
 256       if(nTrack->GetTPCNcls() < 50)continue;
 257       if(pTrack->Chi2perNDF() > 4)continue;
 258       if(nTrack->Chi2perNDF() > 4)continue;
 259
 260       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 261       if(!pTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 262       if(TMath::Abs(dca[1]) > 2) continue;
 263       if(!nTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 264       if(TMath::Abs(dca[1]) > 2) continue;
 265
 266       V0Index[nGoodV0s] = iV0;
 267       V0TrackID[2*nGoodV0s] = pTrack->GetID();
 268       V0TrackID[2*nGoodV0s+1] = nTrack->GetID();
 269
 270       nGoodV0s++;
 271       if(nGoodV0s > 2) break;
 272   }//V0s loop
 273
 274   //Track loop
 275   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 276     AliAODTrack *trk = aod->GetTrack(itr);
 277     if( !trk ) continue;
 278     if(!(trk->TestFilterBit(1<<0))) continue;
 279
 280       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 281       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 282       if(trk->GetTPCNcls() < 50)continue;
 283       if(trk->Chi2perNDF() > 4)continue;
 284       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 285       if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 286       if(TMath::Abs(dca[1]) > 2) continue;
 287
 288       TrackID[nGoodTracks] = trk->GetID();
 289       nGoodTracks++;
 290
 291       if(nGoodTracks > 4) break;
 292   }//Track loop
 293
 294   if(nGoodV0s == 2 && nGoodTracks == 4){
 295         SortArray(TrackID);
 296         SortArray(V0TrackID);
 297         for(Int_t i=0; i<4; i++) if (TrackID[i] != V0TrackID[i]) return;
 298
 299         AliAODv0 *v00 = aod->GetV0(V0Index[0]);
 300         AliAODv0 *v01 = aod->GetV0(V0Index[1]);
 301         fHistK0sMass->Fill(v00->MassK0Short(),v01->MassK0Short());
 302
 303   }
 304
 305   PostData(2, fListHist);
 306
 307 }
 308
 309 //_____________________________________________________________________________
 310 void AliAnalysisTaskUpcK0sK0s::RunAODtree()
 311 {
 312   //input event
 313   AliAODEvent *aod = (AliAODEvent*) InputEvent();
 314   if(!aod) return;
 315
 316   //input data
 317   fDataFilnam = aod->GetHeader()->GetESDFileName();
 318   fEvtNum = aod->GetHeader()->GetEventNumberESDFile();
 319   fRunNum = aod->GetRunNumber();
 320
 321
 322   //Trigger
 323   //Trigger
 324   TString trigger = aod->GetFiredTriggerClasses();
 325
 326   fTrigger[0]  = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
 327   fTrigger[1]  = trigger.Contains("CCUP2-B"); // Double gap
 328   fTrigger[2]  = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
 329
 330   Bool_t isTriggered = kFALSE;
 331   for(Int_t i=0; i<ntrg; i++) {
 332     if( fTrigger[i] ) isTriggered = kTRUE;
 333   }
 334   if(!isTriggered ) return;
 335
 336   //trigger inputs
 337   fL0inputs = aod->GetHeader()->GetL0TriggerInputs();
 338   fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
 339
 340   //Event identification
 341   fPerNum = aod ->GetPeriodNumber();
 342   fOrbNum = aod ->GetOrbitNumber();
 343   fBCrossNum = aod ->GetBunchCrossNumber();
 344
 345   //primary vertex
 346   AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
 347   fVtxContrib = fAODVertex->GetNContributors();
 348
 349   //Tracklets
 350   fNtracklets = aod->GetTracklets()->GetNumberOfTracklets();
 351
 352   //VZERO, ZDC
 353   AliAODVZERO *fV0data = aod ->GetVZEROData();
 354   AliAODZDC *fZDCdata = aod->GetZDCData();
 355
 356   fV0Adecision = fV0data->GetV0ADecision();
 357   fV0Cdecision = fV0data->GetV0CDecision();
 358   fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
 359   fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
 360
 361   Int_t nGoodV0s=0;
 362   Int_t V0Index[3] = {-1,-1,-1};
 363   Int_t V0TrackID[6] = {-1,-1,-1,-1,-1,-1};
 364
 365   Int_t nGoodTracks=0;
 366   Int_t TrackID[5] = {-1,-1,-1,-1,-1};
 367
 368   //V0s loop
 369   for(Int_t iV0=0; iV0<aod ->GetNumberOfV0s(); iV0++) {
 370     AliAODv0 *v0 = aod->GetV0(iV0);
 371     if( !v0 ) continue;
 372     Bool_t lOnFlyStatus = v0->GetOnFlyStatus();
 373     if (lOnFlyStatus) continue;
 374
 375     AliAODTrack *pTrack=(AliAODTrack *)v0->GetDaughter(0); //0->Positive Daughter
 376     AliAODTrack *nTrack=(AliAODTrack *)v0->GetDaughter(1); //1->Negative Daughter
 377     if (!pTrack || !nTrack) continue;
 378
 379     if ( pTrack->Charge() == nTrack->Charge())continue;
 380
 381       if(!(pTrack->TestFilterBit(1<<0))) continue;
 382       if(!(nTrack->TestFilterBit(1<<0))) continue;
 383       if(!(pTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 384       if(!(nTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 385       if(!(pTrack->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 386       if(!(nTrack->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 387       if(pTrack->GetTPCNcls() < 50)continue;
 388       if(nTrack->GetTPCNcls() < 50)continue;
 389       if(pTrack->Chi2perNDF() > 4)continue;
 390       if(nTrack->Chi2perNDF() > 4)continue;
 391
 392       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 393       if(!pTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 394       if(TMath::Abs(dca[1]) > 2) continue;
 395       if(!nTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 396       if(TMath::Abs(dca[1]) > 2) continue;
 397
 398       V0Index[nGoodV0s] = iV0;
 399       V0TrackID[2*nGoodV0s] = pTrack->GetID();
 400       V0TrackID[2*nGoodV0s+1] = nTrack->GetID();
 401       nGoodV0s++;
 402       if(nGoodV0s > 2) break;
 403   }//V0s loop
 404
 405   //Track loop
 406   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 407     AliAODTrack *trk = aod->GetTrack(itr);
 408     if( !trk ) continue;
 409     if(!(trk->TestFilterBit(1<<0))) continue;
 410
 411       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 412       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 413       if(trk->GetTPCNcls() < 50)continue;
 414       if(trk->Chi2perNDF() > 4)continue;
 415       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 416       if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 417       if(TMath::Abs(dca[1]) > 2) continue;
 418
 419       TrackID[nGoodTracks] = trk->GetID();
 420       nGoodTracks++;
 421
 422       if(nGoodTracks > 4) break;
 423   }//Track loop
 424
 425   if(nGoodV0s == 2 && nGoodTracks == 4){
 426         //SortArray(TrackID);
 427         //SortArray(V0TrackID);
 428         //for{Int_t i=0; i<4; i++} if (TrackID[i] != V0TrackID[i]) return;
 429         for(Int_t i=0; i<2; i++){
 430                 AliAODv0 *v0 = aod->GetV0(V0Index[i]);
 431                 AliAODTrack *pTrack=(AliAODTrack *)v0->GetDaughter(0); //0->Positive Daughter
 432                 AliAODTrack *nTrack=(AliAODTrack *)v0->GetDaughter(1); //1->Negative Daughter
 433
 434                 new((*fK0sAODv0s)[i]) AliAODv0(*v0);
 435                 new((*fK0sAODTracks)[2*i]) AliAODTrack(*pTrack);
 436                 new((*fK0sAODTracks)[2*i+1]) AliAODTrack(*nTrack);
 437                 }
 438   fK0sTree ->Fill();
 439
 440   }
 441
 442   PostData(1, fK0sTree);
 443
 444 }//RunAOD
 445
 446 //_____________________________________________________________________________
 447 void AliAnalysisTaskUpcK0sK0s::SortArray(Int_t *dArray) {
 448     for (Int_t i = 3; i > 0; --i) {
 449         for (Int_t j = 0; j < i; ++j) {
 450             if (dArray[j] > dArray[j+1]) {
 451                 Int_t dTemp = dArray[j];
 452                 dArray[j] = dArray[j+1];
 453                 dArray[j+1] = dTemp;
 454             }
 455         }
 456     }
 457 }
 458
 459 //_____________________________________________________________________________
 460 void AliAnalysisTaskUpcK0sK0s::RunESD()
 461 {
 462
 463 }//RunESD
 464
 465 //_____________________________________________________________________________
 466 void AliAnalysisTaskUpcK0sK0s::Terminate(Option_t *)
 467 {
 468
 469   cout<<"Analysis complete."<<endl;
 470 }//Terminate
 471