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