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 = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
 277     if(!trk) AliFatal("Not a standard AOD");
 278     if( !trk ) continue;
 279     if(!(trk->TestFilterBit(1<<0))) continue;
 280
 281       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 282       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 283       if(trk->GetTPCNcls() < 50)continue;
 284       if(trk->Chi2perNDF() > 4)continue;
 285       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 286       if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 287       if(TMath::Abs(dca[1]) > 2) continue;
 288
 289       TrackID[nGoodTracks] = trk->GetID();
 290       nGoodTracks++;
 291
 292       if(nGoodTracks > 4) break;
 293   }//Track loop
 294
 295   if(nGoodV0s == 2 && nGoodTracks == 4){
 296         SortArray(TrackID);
 297         SortArray(V0TrackID);
 298         for(Int_t i=0; i<4; i++) if (TrackID[i] != V0TrackID[i]) return;
 299
 300         AliAODv0 *v00 = aod->GetV0(V0Index[0]);
 301         AliAODv0 *v01 = aod->GetV0(V0Index[1]);
 302         fHistK0sMass->Fill(v00->MassK0Short(),v01->MassK0Short());
 303
 304   }
 305
 306   PostData(2, fListHist);
 307
 308 }
 309
 310 //_____________________________________________________________________________
 311 void AliAnalysisTaskUpcK0sK0s::RunAODtree()
 312 {
 313   //input event
 314   AliAODEvent *aod = (AliAODEvent*) InputEvent();
 315   if(!aod) return;
 316
 317   //input data
 318   AliAODHeader * header = dynamic_cast<AliAODHeader*>(aod->GetHeader());
 319 if(!header) AliFatal("Not a standard AOD");
 320
 321   fDataFilnam = header->GetESDFileName();
 322   fEvtNum = header->GetEventNumberESDFile();
 323   fRunNum = aod->GetRunNumber();
 324
 325
 326   //Trigger
 327   //Trigger
 328   TString trigger = aod->GetFiredTriggerClasses();
 329
 330   fTrigger[0]  = trigger.Contains("CCUP4-B"); // Central UPC Pb-Pb 2011
 331   fTrigger[1]  = trigger.Contains("CCUP2-B"); // Double gap
 332   fTrigger[2]  = trigger.Contains("CCUP7-B"); // Central UPC p-Pb 2013
 333
 334   Bool_t isTriggered = kFALSE;
 335   for(Int_t i=0; i<ntrg; i++) {
 336     if( fTrigger[i] ) isTriggered = kTRUE;
 337   }
 338   if(!isTriggered ) return;
 339
 340   //trigger inputs
 341   fL0inputs = ((AliVAODHeader*)aod->GetHeader())->GetL0TriggerInputs();
 342   fL1inputs = ((AliVAODHeader*)aod->GetHeader())->GetL1TriggerInputs();
 343
 344   //Event identification
 345   fPerNum = aod ->GetPeriodNumber();
 346   fOrbNum = aod ->GetOrbitNumber();
 347   fBCrossNum = aod ->GetBunchCrossNumber();
 348
 349   //primary vertex
 350   AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
 351   fVtxContrib = fAODVertex->GetNContributors();
 352
 353   //Tracklets
 354   fNtracklets = aod->GetTracklets()->GetNumberOfTracklets();
 355
 356   //VZERO, ZDC
 357   AliAODVZERO *fV0data = aod ->GetVZEROData();
 358   AliAODZDC *fZDCdata = aod->GetZDCData();
 359
 360   fV0Adecision = fV0data->GetV0ADecision();
 361   fV0Cdecision = fV0data->GetV0CDecision();
 362   fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
 363   fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
 364
 365   Int_t nGoodV0s=0;
 366   Int_t V0Index[3] = {-1,-1,-1};
 367   Int_t V0TrackID[6] = {-1,-1,-1,-1,-1,-1};
 368
 369   Int_t nGoodTracks=0;
 370   Int_t TrackID[5] = {-1,-1,-1,-1,-1};
 371
 372   //V0s loop
 373   for(Int_t iV0=0; iV0<aod ->GetNumberOfV0s(); iV0++) {
 374     AliAODv0 *v0 = aod->GetV0(iV0);
 375     if( !v0 ) continue;
 376     Bool_t lOnFlyStatus = v0->GetOnFlyStatus();
 377     if (lOnFlyStatus) continue;
 378
 379     AliAODTrack *pTrack=(AliAODTrack *)v0->GetDaughter(0); //0->Positive Daughter
 380     AliAODTrack *nTrack=(AliAODTrack *)v0->GetDaughter(1); //1->Negative Daughter
 381     if (!pTrack || !nTrack) continue;
 382
 383     if ( pTrack->Charge() == nTrack->Charge())continue;
 384
 385       if(!(pTrack->TestFilterBit(1<<0))) continue;
 386       if(!(nTrack->TestFilterBit(1<<0))) continue;
 387       if(!(pTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 388       if(!(nTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 389       if(!(pTrack->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 390       if(!(nTrack->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 391       if(pTrack->GetTPCNcls() < 50)continue;
 392       if(nTrack->GetTPCNcls() < 50)continue;
 393       if(pTrack->Chi2perNDF() > 4)continue;
 394       if(nTrack->Chi2perNDF() > 4)continue;
 395
 396       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 397       if(!pTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 398       if(TMath::Abs(dca[1]) > 2) continue;
 399       if(!nTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 400       if(TMath::Abs(dca[1]) > 2) continue;
 401
 402       V0Index[nGoodV0s] = iV0;
 403       V0TrackID[2*nGoodV0s] = pTrack->GetID();
 404       V0TrackID[2*nGoodV0s+1] = nTrack->GetID();
 405       nGoodV0s++;
 406       if(nGoodV0s > 2) break;
 407   }//V0s loop
 408
 409   //Track loop
 410   for(Int_t itr=0; itr<aod ->GetNumberOfTracks(); itr++) {
 411     AliAODTrack *trk = dynamic_cast<AliAODTrack*>(aod->GetTrack(itr));
 412     if(!trk) AliFatal("Not a standard AOD");
 413     if( !trk ) continue;
 414     if(!(trk->TestFilterBit(1<<0))) continue;
 415
 416       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 417       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 418       if(trk->GetTPCNcls() < 50)continue;
 419       if(trk->Chi2perNDF() > 4)continue;
 420       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 421       if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 422       if(TMath::Abs(dca[1]) > 2) continue;
 423
 424       TrackID[nGoodTracks] = trk->GetID();
 425       nGoodTracks++;
 426
 427       if(nGoodTracks > 4) break;
 428   }//Track loop
 429
 430   if(nGoodV0s == 2 && nGoodTracks == 4){
 431         //SortArray(TrackID);
 432         //SortArray(V0TrackID);
 433         //for{Int_t i=0; i<4; i++} if (TrackID[i] != V0TrackID[i]) return;
 434         for(Int_t i=0; i<2; i++){
 435                 AliAODv0 *v0 = aod->GetV0(V0Index[i]);
 436                 AliAODTrack *pTrack=(AliAODTrack *)v0->GetDaughter(0); //0->Positive Daughter
 437                 AliAODTrack *nTrack=(AliAODTrack *)v0->GetDaughter(1); //1->Negative Daughter
 438
 439                 new((*fK0sAODv0s)[i]) AliAODv0(*v0);
 440                 new((*fK0sAODTracks)[2*i]) AliAODTrack(*pTrack);
 441                 new((*fK0sAODTracks)[2*i+1]) AliAODTrack(*nTrack);
 442                 }
 443   fK0sTree ->Fill();
 444
 445   }
 446
 447   PostData(1, fK0sTree);
 448
 449 }//RunAOD
 450
 451 //_____________________________________________________________________________
 452 void AliAnalysisTaskUpcK0sK0s::SortArray(Int_t *dArray) {
 453     for (Int_t i = 3; i > 0; --i) {
 454         for (Int_t j = 0; j < i; ++j) {
 455             if (dArray[j] > dArray[j+1]) {
 456                 Int_t dTemp = dArray[j];
 457                 dArray[j] = dArray[j+1];
 458                 dArray[j+1] = dTemp;
 459             }
 460         }
 461     }
 462 }
 463
 464 //_____________________________________________________________________________
 465 void AliAnalysisTaskUpcK0sK0s::RunESD()
 466 {
 467
 468 }//RunESD
 469
 470 //_____________________________________________________________________________
 471 void AliAnalysisTaskUpcK0sK0s::Terminate(Option_t *)
 472 {
 473
 474   cout<<"Analysis complete."<<endl;
 475 }//Terminate
 476