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 "TObjString.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   //input file
 129   fDataFilnam = new TObjString();
 130   fDataFilnam->SetString("");
 131
 132     //vertices
 133   fK0sAODv0s = new TClonesArray("AliAODv0", 1000);
 134   fK0sESDv0s = new TClonesArray("AliESDv0", 1000);
 135
 136     //tracks
 137   fK0sAODTracks = new TClonesArray("AliAODTrack", 1000);
 138   fK0sESDTracks = new TClonesArray("AliESDtrack", 1000);
 139
 140   //output tree with K0s candidate events
 141   fK0sTree = new TTree("fK0sTree", "fK0sTree");
 142   fK0sTree ->Branch("fRunNum", &fRunNum, "fRunNum/I");
 143   fK0sTree ->Branch("fPerNum", &fPerNum, "fPerNum/i");
 144   fK0sTree ->Branch("fOrbNum", &fOrbNum, "fOrbNum/i");
 145
 146   fK0sTree ->Branch("fBCrossNum", &fBCrossNum, "fBCrossNum/s");
 147   fK0sTree ->Branch("fTrigger", &fTrigger[0], Form("fTrigger[%i]/O", ntrg));
 148   fK0sTree ->Branch("fL0inputs", &fL0inputs, "fL0inputs/i");
 149   fK0sTree ->Branch("fL1inputs", &fL1inputs, "fL1inputs/i");
 150   fK0sTree ->Branch("fNtracklets", &fNtracklets, "fNtracklets/s");
 151   fK0sTree ->Branch("fVtxContrib", &fVtxContrib, "fVtxContrib/I");
 152   fK0sTree ->Branch("fZDCAenergy", &fZDCAenergy, "fZDCAenergy/D");
 153   fK0sTree ->Branch("fZDCCenergy", &fZDCCenergy, "fZDCCenergy/D");
 154   fK0sTree ->Branch("fV0Adecision", &fV0Adecision, "fV0Adecision/I");
 155   fK0sTree ->Branch("fV0Cdecision", &fV0Cdecision, "fV0Cdecision/I");
 156   fK0sTree ->Branch("fDataFilnam", &fDataFilnam);
 157   fK0sTree ->Branch("fRecoPass", &fRecoPass, "fRecoPass/S");
 158   fK0sTree ->Branch("fEvtNum", &fEvtNum, "fEvtNum/L");
 159   if( fType == 0 ) {
 160     fK0sTree ->Branch("fK0sESDv0s", &fK0sESDv0s);
 161     fK0sTree ->Branch("fK0sESDTracks", &fK0sESDTracks);
 162   }
 163   if( fType == 1 ) {
 164     fK0sTree ->Branch("fK0sAODv0s", &fK0sAODv0s);
 165     fK0sTree ->Branch("fK0sAODTracks", &fK0sAODTracks);
 166   }
 167
 168
 169   fListHist = new TList();
 170   fListHist ->SetOwner();
 171
 172   fHistTriggersPerRun = new TH1D("fHistTriggersPerRun", "fHistTriggersPerRun", 3000, 167000.5, 170000.5);
 173   fListHist->Add(fHistTriggersPerRun);
 174
 175   fHistK0sMass = new TH2D("fHistK0sMass","fHistK0sMass",20,0.45,0.55,20,0.45,0.55);
 176   fListHist->Add(fHistK0sMass);
 177
 178   PostData(1, fK0sTree);
 179   PostData(2, fListHist);
 180
 181 }//UserCreateOutputObjects
 182
 183 //_____________________________________________________________________________
 184 void AliAnalysisTaskUpcK0sK0s::UserExec(Option_t *)
 185 {
 186
 187   //cout<<"#################### Next event ##################"<<endl;
 188
 189   if( fType == 0 ) RunESD();
 190   if( fType == 1 ){
 191         if(fRunHist) RunAODhist();
 192         if(fRunTree) RunAODtree();
 193         }
 194
 195 }//UserExec
 196 //_____________________________________________________________________________
 197 void AliAnalysisTaskUpcK0sK0s::RunAODhist()
 198 {
 199
 200   //input event
 201   AliAODEvent *aod = (AliAODEvent*) InputEvent();
 202   if(!aod) return;
 203
 204
 205   //Trigger
 206   TString trigger = aod->GetFiredTriggerClasses();
 207
 208   if( !trigger.Contains("CCUP4-B") ) return;
 209
 210   fRunNum = aod ->GetRunNumber();
 211   fHistTriggersPerRun->Fill(fRunNum);
 212
 213
 214   //primary vertex
 215   AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
 216   fVtxContrib = fAODVertex->GetNContributors();
 217   if(fVtxContrib < 2) return;
 218
 219
 220   //VZERO, ZDC
 221   AliAODVZERO *fV0data = aod ->GetVZEROData();
 222
 223   fV0Adecision = fV0data->GetV0ADecision();
 224   fV0Cdecision = fV0data->GetV0CDecision();
 225   if(fV0Adecision != AliAODVZERO::kV0Empty || fV0Cdecision != AliAODVZERO::kV0Empty) return;
 226
 227   AliAODZDC *fZDCdata = aod->GetZDCData();
 228
 229   fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
 230   fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
 231   if( fZDCAenergy > 6000 || fZDCCenergy > 6000) return;
 232
 233   Int_t nGoodV0s=0;
 234   Int_t V0Index[3] = {-1,-1,-1};
 235   Int_t V0TrackID[6] = {-1,-1,-1,-1,-1,-1};
 236
 237   Int_t nGoodTracks=0;
 238   Int_t TrackID[5] = {-1,-1,-1,-1,-1};
 239
 240
 241   //V0s loop
 242   for(Int_t iV0=0; iV0<aod ->GetNumberOfV0s(); iV0++) {
 243     AliAODv0 *v0 = aod->GetV0(iV0);
 244     if( !v0 ) continue;
 245     Bool_t lOnFlyStatus = v0->GetOnFlyStatus();
 246     if (lOnFlyStatus) continue;
 247
 248     AliAODTrack *pTrack=(AliAODTrack *)v0->GetDaughter(0); //0->Positive Daughter
 249     AliAODTrack *nTrack=(AliAODTrack *)v0->GetDaughter(1); //1->Negative Daughter
 250     if (!pTrack || !nTrack) continue;
 251
 252     if ( pTrack->Charge() == nTrack->Charge())continue;
 253
 254       if(!(pTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 255       if(!(nTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 256       if(pTrack->GetTPCNcls() < 50)continue;
 257       if(nTrack->GetTPCNcls() < 50)continue;
 258       if(pTrack->Chi2perNDF() > 4)continue;
 259       if(nTrack->Chi2perNDF() > 4)continue;
 260
 261       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 262       if(!pTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 263       if(TMath::Abs(dca[1]) > 2) continue;
 264       if(!nTrack->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 265       if(TMath::Abs(dca[1]) > 2) continue;
 266
 267       V0Index[nGoodV0s] = iV0;
 268       V0TrackID[2*nGoodV0s] = pTrack->GetID();
 269       V0TrackID[2*nGoodV0s+1] = nTrack->GetID();
 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
 279       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 280       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 281       if(trk->GetTPCNcls() < 50)continue;
 282       if(trk->Chi2perNDF() > 4)continue;
 283       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 284       if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 285       if(TMath::Abs(dca[1]) > 2) continue;
 286
 287       TrackID[nGoodTracks] = trk->GetID();
 288       nGoodTracks++;
 289
 290       if(nGoodTracks > 4) break;
 291   }//Track loop
 292
 293   if(nGoodV0s == 2 && nGoodTracks == 4){
 294         //SortArray(TrackID);
 295         //SortArray(V0TrackID);
 296         //for(Int_t i=0; i<4; i++) if (TrackID[i] != V0TrackID[i]) return;
 297         AliAODv0 *v00 = aod->GetV0(V0Index[0]);
 298         AliAODv0 *v01 = aod->GetV0(V0Index[1]);
 299         fHistK0sMass->Fill(v00->MassK0Short(),v01->MassK0Short());
 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   const char *filnam = ((TTree*) GetInputData(0))->GetCurrentFile()->GetName();
 315   fDataFilnam->Clear();
 316   fDataFilnam->SetString(filnam);
 317   fEvtNum = ((TTree*) GetInputData(0))->GetTree()->GetReadEntry();
 318   fRunNum = aod ->GetRunNumber();
 319
 320   //Trigger
 321   TString trigger = aod->GetFiredTriggerClasses();
 322
 323   if( !trigger.Contains("CCUP4-B") ) return;
 324
 325   //trigger inputs
 326   fL0inputs = aod->GetHeader()->GetL0TriggerInputs();
 327   fL1inputs = aod->GetHeader()->GetL1TriggerInputs();
 328
 329   //Event identification
 330   fPerNum = aod ->GetPeriodNumber();
 331   fOrbNum = aod ->GetOrbitNumber();
 332   fBCrossNum = aod ->GetBunchCrossNumber();
 333
 334   //primary vertex
 335   AliAODVertex *fAODVertex = aod->GetPrimaryVertex();
 336   fVtxContrib = fAODVertex->GetNContributors();
 337
 338   //Tracklets
 339   fNtracklets = aod->GetTracklets()->GetNumberOfTracklets();
 340
 341   //VZERO, ZDC
 342   AliAODVZERO *fV0data = aod ->GetVZEROData();
 343   AliAODZDC *fZDCdata = aod->GetZDCData();
 344
 345   fV0Adecision = fV0data->GetV0ADecision();
 346   fV0Cdecision = fV0data->GetV0CDecision();
 347   fZDCAenergy = fZDCdata->GetZNATowerEnergy()[0];
 348   fZDCCenergy = fZDCdata->GetZNCTowerEnergy()[0];
 349
 350   Int_t nGoodV0s=0;
 351   Int_t V0Index[3] = {-1,-1,-1};
 352   Int_t V0TrackID[6] = {-1,-1,-1,-1,-1,-1};
 353
 354   Int_t nGoodTracks=0;
 355   Int_t TrackID[5] = {-1,-1,-1,-1,-1};
 356
 357   //V0s loop
 358   for(Int_t iV0=0; iV0<aod ->GetNumberOfV0s(); iV0++) {
 359     AliAODv0 *v0 = aod->GetV0(iV0);
 360     if( !v0 ) continue;
 361     Bool_t lOnFlyStatus = v0->GetOnFlyStatus();
 362     if (lOnFlyStatus) continue;
 363
 364     AliAODTrack *pTrack=(AliAODTrack *)v0->GetDaughter(0); //0->Positive Daughter
 365     AliAODTrack *nTrack=(AliAODTrack *)v0->GetDaughter(1); //1->Negative Daughter
 366     if (!pTrack || !nTrack) continue;
 367
 368     if ( pTrack->Charge() == nTrack->Charge())continue;
 369
 370       if(!(pTrack->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 371       if(!(nTrack->GetStatus() & AliESDtrack::kTPCrefit) ) 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
 395       if(!(trk->GetStatus() & AliESDtrack::kTPCrefit) ) continue;
 396       if(!(trk->GetStatus() & AliESDtrack::kITSrefit) ) continue;
 397       if(trk->GetTPCNcls() < 50)continue;
 398       if(trk->Chi2perNDF() > 4)continue;
 399       Double_t dca[2] = {0.0,0.0}, cov[3] = {0.0,0.0,0.0};
 400       if(!trk->PropagateToDCA(fAODVertex,aod->GetMagneticField(),300.,dca,cov)) continue;
 401       if(TMath::Abs(dca[1]) > 2) continue;
 402
 403       TrackID[nGoodTracks] = trk->GetID();
 404       nGoodTracks++;
 405
 406       if(nGoodTracks > 4) break;
 407   }//Track loop
 408
 409   if(nGoodV0s == 2 && nGoodTracks == 4){
 410         //SortArray(TrackID);
 411         //SortArray(V0TrackID);
 412         //for{Int_t i=0; i<4; i++} if (TrackID[i] != V0TrackID[i]) return;
 413         for(Int_t i=0; i<2; i++){
 414                 AliAODv0 *v0 = aod->GetV0(V0Index[i]);
 415                 AliAODTrack *pTrack=(AliAODTrack *)v0->GetDaughter(0); //0->Positive Daughter
 416                 AliAODTrack *nTrack=(AliAODTrack *)v0->GetDaughter(1); //1->Negative Daughter
 417
 418                 new((*fK0sAODv0s)[i]) AliAODv0(*v0);
 419                 new((*fK0sAODTracks)[2*i]) AliAODTrack(*pTrack);
 420                 new((*fK0sAODTracks)[2*i+1]) AliAODTrack(*nTrack);
 421                 }
 422   fK0sTree ->Fill();
 423
 424   }
 425
 426   PostData(1, fK0sTree);
 427
 428 }//RunAOD
 429
 430 //_____________________________________________________________________________
 431 void AliAnalysisTaskUpcK0sK0s::SortArray(Double_t *dArray) {
 432     for (Int_t i = 3; i > 0; --i) {
 433         for (Int_t j = 0; j < i; ++j) {
 434             if (dArray[j] > dArray[j+1]) {
 435                 Double_t dTemp = dArray[j];
 436                 dArray[j] = dArray[j+1];
 437                 dArray[j+1] = dTemp;
 438             }
 439         }
 440     }
 441 }
 442
 443 //_____________________________________________________________________________
 444 void AliAnalysisTaskUpcK0sK0s::RunESD()
 445 {
 446
 447 }//RunESD
 448
 449 //_____________________________________________________________________________
 450 void AliAnalysisTaskUpcK0sK0s::Terminate(Option_t *)
 451 {
 452
 453   cout<<"Analysis complete."<<endl;
 454 }//Terminate
 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
 484