PWG4/GammaConv/AliAnalysisTaskGCPartToPWG4Part.cxx

   1 #include <iostream>
   2 #include "TChain.h"
   3 #include "TTree.h"
   4 #include "TH1F.h"
   5 #include "TCanvas.h"
   6 #include "TString.h"
   7
   8 #include "AliAnalysisTask.h"
   9 #include "AliAnalysisManager.h"
  10 #include "AliAnalysisTaskGCPartToPWG4Part.h"
  11
  12 #include "AliESDEvent.h"
  13 #include "AliESDCaloCluster.h"
  14 #include "AliESDInputHandler.h"
  15
  16 #include "AliAODPWG4ParticleCorrelation.h"
  17 #include "AliAODEvent.h"
  18 #include "AliAODHandler.h"
  19 #include "AliAODCaloCluster.h"
  20 #include "AliGammaConversionAODObject.h"
  21 #include "AliAODConversionParticle.h"
  22 #include "AliAODJet.h"
  23
  24 #include "AliAODInputHandler.h"
  25
  26 #include "AliAODMCParticle.h"
  27
  28 #include "AliMCAnalysisUtils.h"
  29
  30
  31 #include "AliAODMCHeader.h"
  32 // Gamma - jet correlation analysis task
  33 // Authors: Svein Lindal
  34
  35
  36 using namespace std;
  37
  38 ClassImp(AliAnalysisTaskGCPartToPWG4Part)
  39
  40 //________________________________________________________________________
  41 AliAnalysisTaskGCPartToPWG4Part::AliAnalysisTaskGCPartToPWG4Part()
  42 : AliAnalysisTaskSE(),
  43   fDeltaAODFileName(""),
  44   fAODBranchName("GammaConv_gamma"),
  45   fAODPWG4Particles(NULL),
  46   fAnaUtils(NULL)
  47 {
  48   // Dummy Constructor
  49 }
  50
  51 //________________________________________________________________________________
  52 AliAnalysisTaskGCPartToPWG4Part::~AliAnalysisTaskGCPartToPWG4Part() {
  53
  54   if(fAODPWG4Particles)
  55     fAODPWG4Particles = NULL;
  56   delete fAODPWG4Particles;
  57
  58   if(fAnaUtils)
  59     delete fAnaUtils;
  60   fAnaUtils = NULL;
  61
  62 }
  63
  64
  65
  66 //________________________________________________________________________
  67 AliAnalysisTaskGCPartToPWG4Part::AliAnalysisTaskGCPartToPWG4Part(const char *name) :
  68   AliAnalysisTaskSE(name),
  69   fDeltaAODFileName(""),
  70   fAODBranchName("GammaConv_gamma"),
  71   fAODPWG4Particles(NULL),
  72   fAnaUtils(NULL)
  73 {
  74   // Constructor
  75   DefineInput(0, TChain::Class());
  76   // Output slot #0 id reserved by the base class for AOD
  77
  78   // Output slot #1 writes into a TH1 container
  79   DefineOutput(1, TList::Class());
  80
  81   fAnaUtils = new AliMCAnalysisUtils();
  82 }
  83
  84
  85
  86 //________________________________________________________________________
  87 void AliAnalysisTaskGCPartToPWG4Part::UserCreateOutputObjects() {
  88   fAODPWG4Particles = new TClonesArray("AliAODPWG4ParticleCorrelation", 0);
  89   fAODPWG4Particles->SetName("ConversionGamma");
  90   AddAODBranch("TClonesArray", &fAODPWG4Particles);
  91
  92 }
  93
  94 //________________________________________________________________________
  95 void AliAnalysisTaskGCPartToPWG4Part::UserExec(Option_t *)
  96 {
  97
  98   //AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
  99
 100   //Clear stuff for new event
 101   CleanUp();
 102
 103   ///Get AOD event
 104   AliAODEvent * aodEvent = GetAODEvent();
 105   if(!aodEvent) {
 106     AliError("No AOD event!!");
 107     return;
 108   }
 109
 110
 111   ProcessConvGamma(aodEvent);
 112
 113
 114   //PostData(1, fOutputList);
 115
 116 }
 117
 118
 119
 120
 121 //___________________________________________________________________________________________
 122 void AliAnalysisTaskGCPartToPWG4Part::ProcessConvGamma( const AliAODEvent * const aodEvent ) {
 123
 124   TClonesArray * tracks = aodEvent->GetTracks();
 125   if(!tracks) {
 126     cout << "No tracks!!!"<<endl;
 127     return;
 128   }
 129
 130
 131   TClonesArray * convGamma = GetConversionGammas(aodEvent);
 132   if(!convGamma) {
 133     AliError(Form("No branch by name %s found in file %s", fAODBranchName.Data(), fDeltaAODFileName.Data()));
 134     return;
 135   }
 136
 137
 138
 139   // AliAODMCHeader * mcHeader = dynamic_cast<AliAODMCHeader*>(aodEvent->FindListObject("mcHeader"));
 140
 141
 142   TClonesArray * arrayMC = dynamic_cast<TClonesArray*>(aodEvent->GetList()->FindObject(AliAODMCParticle::StdBranchName()));
 143
 144 for (Int_t iPhot = 0; iPhot < convGamma->GetEntriesFast(); iPhot++) {
 145     //AliAODPWG4Particle * photon = dynamic_cast<AliAODPWG4Particle*>(convGamma->At(iPhot));
 146
 147
 148     //if(!photon) {
 149       // AliGammaConversionAODObject * aodO = dynamic_cast<AliGammaConversionAODObject*>(convGamma->At(iPhot));
 150       // if (!aodO) {
 151
 152       //        AliError(Form("ERROR: Could not receive ga %d\n", iPhot));
 153       //        continue;
 154       // }
 155
 156     AliAODConversionParticle * convParticle = dynamic_cast<AliAODConversionParticle*>(convGamma->At(iPhot));
 157     if (!convParticle) {
 158
 159       AliError(Form("ERROR: Could not receive ga %d\n", iPhot));
 160       continue;
 161     }
 162
 163
 164     AliAODPWG4ParticleCorrelation * photon = AddToAOD(convParticle, fAODPWG4Particles, "ConvGamma");
 165       //Int_t tag = CheckTag(photon, tracks, arrayMC, mcHeader);
 166       //if(tag > 0 ) photon->SetTag(tag);
 167
 168   }
 169
 170
 171 }
 172
 173 /////____________________________________________________________________________________
 174 // void AliAnalysisTaskGCPartToPWG4Part::FillMCHistograms(AliAODPWG4ParticleCorrelation * recParticle, AliAODMCParticle * mcParticle) {
 175
 176
 177
 178
 179 // }
 180 //////_________________________________________________________________________________________
 181 Int_t AliAnalysisTaskGCPartToPWG4Part::CheckTag(AliAODPWG4ParticleCorrelation * particle, TClonesArray * tracks, TClonesArray * arrayMC, AliAODMCHeader * mcHeader) {
 182
 183
 184
 185
 186
 187   for (int imc = 0; imc < arrayMC->GetEntriesFast(); imc++) {
 188   //for (int imc = 0; imc < 20; imc++) {
 189
 190
 191
 192     AliAODMCParticle * mParticle = dynamic_cast<AliAODMCParticle*>(arrayMC->At(imc));
 193     //    cout << mParticle->GetPdgCode() << " " <<mParticle->GetStatus() << endl;
 194
 195
 196     //if( mParticle->GetPdgCode() == 34) cout<< "BALLE"<<endl;// && mParticle->GetMother() < 10) {
 197
 198     if( mParticle->GetPdgCode() == 22 ) { //&& mParticle->GetMother() < 10) {
 199       AliAODMCParticle * mother = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mParticle->GetMother()));
 200       //if (TMath::Abs(mother->GetPdgCode()) < 22 && TMath::Abs(mother->GetPdgCode()) != 11) cout  <<"KKKKKKKKK"<< mother->GetPdgCode() <<endl;
 201
 202       for (int id = 0; id < mParticle->GetNDaughters(); id++) {
 203         int idaughter = mParticle->GetDaughter(id);
 204         if(idaughter>0 && idaughter < arrayMC->GetEntriesFast()) {
 205           AliAODMCParticle * dmParticle = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mParticle->GetDaughter(id)));
 206           //cout  << dmParticle->GetPdgCode() << endl;
 207         }
 208       }
 209     }
 210   }
 211
 212
 213
 214
 215
 216
 217
 218
 219
 220
 221
 222
 223   Int_t tag = 0;
 224
 225   Int_t l1 = particle->GetTrackLabel(0);
 226   Int_t l2 = particle->GetTrackLabel(1);
 227
 228   AliAODTrack * track1 = NULL;
 229   AliAODTrack * track2 = NULL;
 230
 231
 232   for(int i = 0; i < tracks->GetEntriesFast(); i++) {
 233
 234     AliAODTrack * track = (AliAODTrack*)tracks->At(i);
 235     if (track->GetID() == l1) {
 236       track1 = track;
 237     } else if (track->GetID() == l2) {
 238       track2 = track;
 239     }
 240
 241     if(track1 && track2) break;
 242   }
 243
 244
 245   if(!track1 || !track2) return tag;
 246   if(track1->GetLabel() < 0 || track2->GetLabel() < 0) {
 247     //cout << "error balla"<< endl;
 248
 249   } else {
 250
 251     AliAODMCParticle * mcPart1 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(track1->GetLabel()));
 252     AliAODMCParticle * mcPart2 = dynamic_cast<AliAODMCParticle*>(arrayMC->At(track2->GetLabel()));
 253
 254     if (mcPart1 && mcPart2) {
 255
 256       if(mcPart1->GetMother() == mcPart2->GetMother()) {
 257
 258         AliAODMCParticle * photon = dynamic_cast<AliAODMCParticle*>(arrayMC->At(mcPart1->GetMother()));
 259         Int_t motherIndex = photon->GetMother();
 260         tag= fAnaUtils->CheckOriginInAOD(&motherIndex, 1, arrayMC);
 261
 262         AliAODMCParticle * mother = dynamic_cast<AliAODMCParticle*>(arrayMC->At(motherIndex));
 263
 264
 265         if (photon->GetPdgCode() == 22 ) {
 266
 267           if (TMath::Abs(mother->GetPdgCode()) < 22 && TMath::Abs(mother->GetPdgCode()) != 11) {
 268
 269             fAnaUtils->SetTagBit(tag, AliMCAnalysisUtils::kMCPhoton);
 270             fAnaUtils->SetTagBit(tag, AliMCAnalysisUtils::kMCPrompt);
 271
 272             cout  <<"KKKKKKKKK "<< mother->GetPdgCode() << " " << mother->GetStatus() << " mi " << motherIndex << " daught: " << mother->GetNDaughters() << endl;
 273             cout << photon->GetStatus() << " " << photon->Pt() << " headerpt:" << mcHeader->GetPtHard() << " " <<mcPart1->GetMother() << endl;
 274
 275           }
 276         } else {
 277           cout << "FAKE "<< photon->PdgCode() << " " << photon->Pt() << " " <<  photon->GetStatus() << endl;
 278         }
 279
 280
 281
 282
 283
 284
 285         Int_t parentId = mother->GetMother();
 286         AliAODMCParticle * gp = dynamic_cast<AliAODMCParticle*>(arrayMC->At(parentId));
 287
 288
 289
 290
 291         if(! mother->IsPrimary()) {
 292           //cout << mother->GetPdgCode() << " " << mother->GetStatus() <<endl;
 293           //      cout << "other one? " << mother->PdgCode() << " " << gp->PdgCode() << " " << gp->GetStatus() << endl;
 294           //    cout << "yeay3  " <<  mother->GetPdgCode() << endl;
 295
 296         } else { //if (mother->IsPhysicalPrimary() ){
 297
 298
 299           if(mother->GetStatus() < 10) {
 300             //cout << "yeay4  " <<  mother->GetPdgCode() << endl;
 301             // cout << "pion? " << gp->PdgCode() << " " << gp->GetStatus() << endl;
 302             //cout << "???" << mother->PdgCode() << " " << mother->GetStatus() << " " << gp->PdgCode() << " " << gp->GetStatus() << endl;
 303           } else {
 304             //cout << "yeay5  " <<  mother->GetPdgCode() << endl;
 305             //cout << "other? " << mother->PdgCode() << " " << mother->GetStatus() << " " << gp->PdgCode() << " " << gp->GetStatus() << endl;
 306           }
 307
 308
 309           // } else {
 310           //   cout << "MI: " << mother->GetPdgCode() << " " << mother->GetStatus() <<endl;
 311         }
 312
 313       }
 314     }
 315   }
 316   //cout << "REturn tag " << tag << endl;
 317   return tag;
 318 }
 319
 320
 321 ///__________________________________________________________________________________
 322 AliAODPWG4ParticleCorrelation * AliAnalysisTaskGCPartToPWG4Part::AddToAOD(AliGammaConversionAODObject * aodO, TClonesArray * branch, TString detector) {
 323   new((*branch)[branch->GetEntriesFast()]) AliAODPWG4ParticleCorrelation(aodO->Px(), aodO->Py(), aodO->Pz(), aodO->E());
 324   AliAODPWG4ParticleCorrelation * photon = dynamic_cast<AliAODPWG4ParticleCorrelation*>(branch->Last());
 325   photon->SetTagged(aodO->IsTagged());
 326   photon->SetTrackLabel(aodO->GetLabel1(), aodO->GetLabel2());
 327   photon->SetDetector(detector);
 328   return photon;
 329 }
 330
 331 ///__________________________________________________________________________________
 332 AliAODPWG4ParticleCorrelation * AliAnalysisTaskGCPartToPWG4Part::AddToAOD(AliAODConversionParticle * aodO, TClonesArray * branch, TString detector) {
 333   new((*branch)[branch->GetEntriesFast()]) AliAODPWG4ParticleCorrelation(aodO->Px(), aodO->Py(), aodO->Pz(), aodO->E());
 334   AliAODPWG4ParticleCorrelation * photon = dynamic_cast<AliAODPWG4ParticleCorrelation*>(branch->Last());
 335   photon->SetTrackLabel(aodO->GetLabel1(), aodO->GetLabel2());
 336   photon->SetDetector(detector);
 337   return photon;
 338 }
 339
 340
 341
 342
 343 //_____________________________________________________________________
 344 void AliAnalysisTaskGCPartToPWG4Part::Terminate(Option_t *) {
 345   // Draw result to the screen
 346   // Called once at the end of the query
 347 }
 348
 349 //_____________________________________________________________________
 350 AliAODEvent * AliAnalysisTaskGCPartToPWG4Part::GetAODEvent() {
 351   //Get the AOD event from whereever it might be
 352   AliAODEvent * aodEvent = dynamic_cast<AliAODEvent*>(InputEvent());
 353   if(!aodEvent) {
 354     aodEvent = AODEvent();
 355   }
 356
 357   return aodEvent;
 358
 359 }
 360
 361 //_____________________________________________________________________
 362 TClonesArray * AliAnalysisTaskGCPartToPWG4Part::GetConversionGammas(const AliAODEvent * aodEvent) {
 363
 364   //Get Conversion gamma branch of AOD. First try standard AOD
 365   TClonesArray * convGamma = dynamic_cast<TClonesArray*>(aodEvent->FindListObject(fAODBranchName.Data()));
 366
 367
 368   //If it's there, send it back
 369   if(convGamma)  return convGamma;
 370
 371
 372   //If AOD not in standard file have to locate it in delta AOD
 373   if( !(fDeltaAODFileName.Length() > 0)  ) return NULL;
 374
 375   AliAODHandler * aodHandler = dynamic_cast<AliAODHandler*>(AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler());
 376   if(aodHandler) {
 377     AliAODExtension * gExt = dynamic_cast<AliAODExtension*>(aodHandler->GetExtensions()->FindObject(fDeltaAODFileName));
 378     if(gExt) {
 379       AliAODEvent * gcEvent = gExt->GetAOD();
 380       return dynamic_cast<TClonesArray*>(gcEvent->FindListObject("GammaConv_gamma"));
 381     }
 382   }
 383   return NULL;
 384 }
 385
 386 //_________________________________________________________________________
 387 void AliAnalysisTaskGCPartToPWG4Part::CleanUp() {
 388   fAODPWG4Particles->Delete();
 389 }
 390
 391