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