PWGCF/FEMTOSCOPY/AliFemto/AliAnalysisTaskFemto.cxx

   1 //------------------------------------------------------
   2 // AliAnalysisTaskFemto - A task for the analysis framework
   3 // from the FEMTOSCOPY analysis of PWG2. Creates the necessary
   4 // connection between the ESD or AOD input and the femtoscopic
   5 // code.
   6 // Author: Adam Kisiel, OSU; Adam.Kisiel@cern.ch
   7 //------------------------------------------------------
   8 #include "TROOT.h"
   9 #include "TChain.h"
  10 #include "TH1.h"
  11 #include "TCanvas.h"
  12 #include "TSystem.h"
  13 #include "TFile.h"
  14 #include "TInterpreter.h"
  15
  16 //#include "AliAnalysisTask.h"
  17 #include "AliAnalysisTaskSE.h"
  18
  19 #include "AliESDEvent.h"
  20
  21 #include "AliFemtoAnalysis.h"
  22 #include "AliAnalysisTaskFemto.h"
  23 #include "AliVHeader.h"
  24 #include "AliGenEventHeader.h"
  25 #include "AliGenHijingEventHeader.h"
  26 #include "AliGenCocktailEventHeader.h"
  27
  28 ClassImp(AliAnalysisTaskFemto)
  29
  30 // Default name for the setup macro of femto analysis
  31 // This function MUST be defined in the separate file !!!
  32 // extern AliFemtoManager *ConfigFemtoAnalysis();
  33
  34 //________________________________________________________________________
  35 AliAnalysisTaskFemto::AliAnalysisTaskFemto(const char *name, const char *aConfigMacro, const char *aConfigParams):
  36     AliAnalysisTaskSE(name), //AliAnalysisTask(name,""),
  37     fESD(0),
  38     fESDpid(0),
  39     fAOD(0),
  40     fAODpidUtil(0),
  41     fStack(0),
  42     fOutputList(0),
  43     fReader(0x0),
  44     fManager(0x0),
  45     fAnalysisType(0),
  46     fConfigMacro(0),
  47     fConfigParams(0)
  48 {
  49   // Constructor.
  50   // Input slot #0 works with an Ntuple
  51   DefineInput(0, TChain::Class());
  52   // Output slot #0 writes into a TH1 container
  53   DefineOutput(0, TList::Class());
  54   fConfigMacro = (char *) malloc(sizeof(char) * strlen(aConfigMacro));
  55   strcpy(fConfigMacro, aConfigMacro);
  56   fConfigParams = (char *) malloc(sizeof(char) * strlen(aConfigParams));
  57   strcpy(fConfigParams, aConfigParams);
  58 }
  59 //________________________________________________________________________
  60 AliAnalysisTaskFemto::AliAnalysisTaskFemto(const char *name, const char *aConfigMacro="ConfigFemtoAnalysis.C"):
  61     AliAnalysisTaskSE(name), //AliAnalysisTask(name,""),
  62     fESD(0),
  63     fESDpid(0),
  64     fAOD(0),
  65     fAODpidUtil(0),
  66     fStack(0),
  67     fOutputList(0),
  68     fReader(0x0),
  69     fManager(0x0),
  70     fAnalysisType(0),
  71     fConfigMacro(0),
  72     fConfigParams(0)
  73 {
  74   // Constructor.
  75   // Input slot #0 works with an Ntuple
  76   DefineInput(0, TChain::Class());
  77   // Output slot #0 writes into a TH1 container
  78   DefineOutput(0, TList::Class());
  79   fConfigMacro = (char *) malloc(sizeof(char) * strlen(aConfigMacro));
  80   strcpy(fConfigMacro, aConfigMacro);
  81   fConfigParams = (char *) malloc(sizeof(char) * 2);
  82   strcpy(fConfigParams, "");
  83 }
  84
  85 AliAnalysisTaskFemto::AliAnalysisTaskFemto(const AliAnalysisTaskFemto& aFemtoTask):
  86     AliAnalysisTaskSE(aFemtoTask), //AliAnalysisTask(aFemtoTask),
  87     fESD(0),
  88     fESDpid(0),
  89     fAOD(0),
  90     fAODpidUtil(0),
  91     fStack(0),
  92     fOutputList(0),
  93     fReader(0x0),
  94     fManager(0x0),
  95     fAnalysisType(0),
  96     fConfigMacro(0),
  97     fConfigParams(0)
  98 {
  99   // copy constructor
 100   fESD = aFemtoTask.fESD;
 101   fESDpid = aFemtoTask.fESDpid;
 102   fAOD = aFemtoTask.fAOD;
 103   fAODpidUtil = aFemtoTask.fAODpidUtil;
 104   fStack = aFemtoTask.fStack;
 105   fOutputList = aFemtoTask.fOutputList;
 106   fReader = aFemtoTask.fReader;
 107   fManager = aFemtoTask.fManager;
 108   fAnalysisType = aFemtoTask.fAnalysisType;
 109   fConfigMacro = (char *) malloc(sizeof(char) * strlen(aFemtoTask.fConfigMacro));
 110   strcpy(fConfigMacro, aFemtoTask.fConfigMacro);
 111   fConfigParams = (char *) malloc(sizeof(char) * strlen(aFemtoTask.fConfigParams));
 112   strcpy(fConfigParams, aFemtoTask.fConfigParams);
 113 }
 114
 115
 116 AliAnalysisTaskFemto& AliAnalysisTaskFemto::operator=(const AliAnalysisTaskFemto& aFemtoTask){
 117   // assignment operator
 118   if (this == &aFemtoTask)
 119     return *this;
 120
 121   fESD = aFemtoTask.fESD;
 122   fESDpid = aFemtoTask.fESDpid;
 123   fAOD = aFemtoTask.fAOD;
 124   fAODpidUtil = aFemtoTask.fAODpidUtil;
 125   fStack = aFemtoTask.fStack;
 126   fOutputList = aFemtoTask.fOutputList;
 127   fReader = aFemtoTask.fReader;
 128   fManager = aFemtoTask.fManager;
 129   fAnalysisType = aFemtoTask.fAnalysisType;
 130   if (fConfigMacro) free(fConfigMacro);
 131   fConfigMacro = (char *) malloc(sizeof(char) * strlen(aFemtoTask.fConfigMacro));
 132   strcpy(fConfigMacro, aFemtoTask.fConfigMacro);
 133   if (fConfigParams) free(fConfigParams);
 134   fConfigParams = (char *) malloc(sizeof(char) * strlen(aFemtoTask.fConfigParams));
 135   strcpy(fConfigParams, aFemtoTask.fConfigParams);
 136
 137   return *this;
 138 }
 139
 140 AliAnalysisTaskFemto::~AliAnalysisTaskFemto()
 141 {
 142   if (fConfigMacro) free(fConfigMacro);
 143   if (fConfigParams) free(fConfigParams);
 144 }
 145
 146
 147 //________________________________________________________________________
 148 void AliAnalysisTaskFemto::ConnectInputData(Option_t *) {
 149   AliInfo(Form("   ConnectInputData %s\n", GetName()));
 150
 151   fESD = 0;
 152   fESDpid = 0;
 153   fAOD = 0;
 154   fAODpidUtil = 0;
 155   fAnalysisType = 0;
 156
 157   TTree* tree = dynamic_cast<TTree*> (GetInputData(0));
 158   if (!tree) {
 159     AliWarning("Could not read chain from input slot 0");
 160     return;
 161   }
 162
 163   AliFemtoEventReaderESDChain *femtoReader = dynamic_cast<AliFemtoEventReaderESDChain *> (fReader);
 164   if ((dynamic_cast<AliFemtoEventReaderESDChain *> (fReader))) {
 165     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 166
 167     if(esdH) {
 168       AliInfo("Selected ESD analysis");
 169       fAnalysisType = 1;
 170
 171 //       if (!esdH) {
 172 //      AliWarning("Could not get ESDInputHandler");
 173 //       }
 174 //       else {
 175         fESD = esdH->GetEvent();
 176         fESDpid = esdH->GetESDpid();
 177         femtoReader->SetESDPid(fESDpid);
 178 //       }
 179     }
 180 }
 181   else if ((dynamic_cast<AliFemtoEventReaderKinematicsChain *> (fReader))) {
 182     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 183
 184     if(esdH) {
 185       AliInfo("Selected ESD analysis");
 186       fAnalysisType = 1;
 187
 188 //       if (!esdH) {
 189 //      AliWarning("Could not get ESDInputHandler");
 190 //       }
 191 //       else {
 192         fESD = esdH->GetEvent();
 193         //fESDpid = esdH->GetESDpid();
 194         //femtoReader->SetESDPid(fESDpid);
 195 //       }
 196     }
 197   }
 198
 199
 200   AliFemtoEventReaderESDChainKine *femtoReaderESDKine = dynamic_cast<AliFemtoEventReaderESDChainKine *> (fReader);
 201   if ((dynamic_cast<AliFemtoEventReaderESDChainKine *> (fReader))) {
 202     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 203
 204     if(esdH) {
 205       AliInfo("Selected ESD analysis");
 206       fAnalysisType = 1;
 207
 208 //       if (!esdH) {
 209 //      AliWarning("Could not get ESDInputHandler");
 210 //       }
 211 //       else {
 212         fESD = esdH->GetEvent();
 213         fESDpid = esdH->GetESDpid();
 214         femtoReaderESDKine->SetESDPid(fESDpid);
 215 //       }
 216     }
 217 }
 218
 219   //  AliFemtoEventReaderKinematicsChain *femtoReaderKine = dynamic_cast<AliFemtoEventReaderKinematicsChain *> (fReader);
 220   if ((dynamic_cast<AliFemtoEventReaderKinematicsChain *> (fReader))) {
 221     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 222
 223     if(esdH) {
 224       AliInfo("Selected ESD analysis");
 225       fAnalysisType = 1;
 226
 227       //       if (!esdH) {
 228       //        AliWarning("Could not get ESDInputHandler");
 229       //       }
 230       //       else {
 231       fESD = esdH->GetEvent();
 232       //fESDpid = esdH->GetESDpid();
 233       //femtoReader->SetESDPid(fESDpid);
 234       //       }
 235     }
 236   }
 237
 238
 239   AliFemtoEventReaderAODChain *femtoReaderAOD = dynamic_cast<AliFemtoEventReaderAODChain *> (fReader);
 240   if (dynamic_cast<AliFemtoEventReaderAODChain *> (fReader)) {
 241     AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 242
 243     if (!aodH) {
 244       TObject *handler = AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler();
 245       AliInfo("Has output handler ");
 246       if( handler && handler->InheritsFrom("AliAODHandler") ) {
 247         AliInfo("Selected AOD analysis");
 248
 249         fAOD = ((AliAODHandler*)handler)->GetAOD();
 250         fAnalysisType = 2;
 251       }
 252       else {
 253         AliWarning("Selected AOD reader but no AOD handler found");
 254       }
 255     }
 256     else {
 257       AliInfo("Selected AOD analysis");
 258       fAnalysisType = 2;
 259
 260       fAOD = aodH->GetEvent();
 261
 262       fAODpidUtil = aodH->GetAODpidUtil();
 263       //      printf("aodH->GetAODpidUtil(): %x",aodH->GetAODpidUtil());
 264       femtoReaderAOD->SetAODpidUtil(fAODpidUtil);
 265     }
 266   }
 267
 268   if ((!fAOD) && (!fESD)) {
 269     AliWarning("Wrong analysis type: Only ESD and AOD types are allowed!");
 270   }
 271 }
 272
 273 //________________________________________________________________________
 274 void AliAnalysisTaskFemto::CreateOutputObjects() {
 275   AliInfo("Creating Femto Analysis objects\n");
 276
 277   gSystem->SetIncludePath("-I$ROOTSYS/include -I./STEERBase/ -I./ESD/ -I./AOD/ -I./ANALYSIS/ -I./ANALYSISalice/ -I./PWG2AOD/AOD -I./PWG2femtoscopy/FEMTOSCOPY/AliFemto -I./PWG2femtoscopyUser/FEMTOSCOPY/AliFemtoUser");
 278   //  char fcm[2000];
 279 //   sprintf(fcm, "%s++", fConfigMacro);
 280 //   gROOT->LoadMacro(fcm);
 281   gROOT->LoadMacro(fConfigMacro);
 282   //  fJetFinder = (AliJetFinder*) gInterpreter->ProcessLine("ConfigJetAnalysis()");
 283   if (!fConfigParams)
 284     SetFemtoManager((AliFemtoManager *) gInterpreter->ProcessLine("ConfigFemtoAnalysis()"));
 285   else
 286     SetFemtoManager((AliFemtoManager *) gInterpreter->ProcessLine(Form("ConfigFemtoAnalysis(%s)", fConfigParams)));
 287
 288   TList *tOL;
 289   fOutputList = fManager->Analysis(0)->GetOutputList();
 290
 291   for (unsigned int ian = 1; ian<fManager->AnalysisCollection()->size(); ian++) {
 292     tOL = fManager->Analysis(ian)->GetOutputList();
 293
 294     TIter nextListCf(tOL);
 295     while (TObject *obj = nextListCf()) {
 296       fOutputList->Add(obj);
 297     }
 298
 299     delete tOL;
 300   }
 301
 302   PostData(0, fOutputList);
 303 }
 304
 305 //________________________________________________________________________
 306 void AliAnalysisTaskFemto::Exec(Option_t *) {
 307   // Task making a femtoscopic analysis.
 308
 309   if (fAnalysisType==1) {
 310     if (!fESD) {
 311       AliWarning("fESD not available");
 312       return;
 313     }
 314
 315     //Get MC data
 316     AliMCEventHandler*    mctruth = (AliMCEventHandler*)
 317       ((AliAnalysisManager::GetAnalysisManager())->GetMCtruthEventHandler());
 318
 319     AliGenHijingEventHeader *hdh = 0;
 320     if(mctruth) {
 321       fStack = mctruth->MCEvent()->Stack();
 322
 323       AliGenCocktailEventHeader *hd = dynamic_cast<AliGenCocktailEventHeader *> (mctruth->MCEvent()->GenEventHeader());
 324
 325       if (hd) {
 326
 327         //      AliInfo ("Got MC cocktail event header %p\n", (void *) hd);
 328         TList *lhd = hd->GetHeaders();
 329         //      AliInfo ("Got list of headers %d\n", lhd->GetEntries());
 330
 331         for (int iterh=0; iterh<lhd->GetEntries(); iterh++)
 332           {
 333             hdh = dynamic_cast<AliGenHijingEventHeader *> (lhd->At(iterh));
 334             //      AliInfo ("HIJING header at %i is %p\n", iterh, (void *) hdh);
 335           }
 336       }
 337     }
 338
 339     // Get ESD
 340     AliESDInputHandler *esdH = dynamic_cast<AliESDInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 341
 342     if (!esdH) {
 343       AliWarning("Could not get ESDInputHandler");
 344       return;
 345     }
 346     else {
 347       fESD = esdH->GetEvent();
 348       fESDpid = esdH->GetESDpid();
 349     }
 350
 351     AliInfo(Form("Tracks in ESD: %d \n",fESD->GetNumberOfTracks()));
 352
 353     if (fESD->GetNumberOfTracks() >= 0) {
 354
 355       if (!fReader) {
 356         AliWarning("No ESD reader for ESD analysis !\n");
 357       }
 358
 359       AliFemtoEventReaderESDChain* fesdc = dynamic_cast<AliFemtoEventReaderESDChain *> (fReader);
 360       if (fesdc)
 361         {
 362           // Process the event with no Kine information
 363           fesdc->SetESDSource(fESD);
 364           fManager->ProcessEvent();
 365         }
 366
 367         AliFemtoEventReaderKinematicsChain* fkinec = dynamic_cast<AliFemtoEventReaderKinematicsChain *> (fReader);
 368       if (fkinec)
 369         {
 370           // Process the event with Kine information only
 371           fkinec->SetStackSource(fStack);
 372           fManager->ProcessEvent();
 373         }
 374
 375
 376       AliFemtoEventReaderESDChainKine* fesdck = dynamic_cast<AliFemtoEventReaderESDChainKine *> (fReader);
 377       if (fesdck)
 378         {
 379           // Process the event with Kine information
 380           fesdck->SetESDSource(fESD);
 381           fesdck->SetStackSource(fStack);
 382
 383           fesdck->SetGenEventHeader(hdh);
 384           fManager->ProcessEvent();
 385         }
 386       AliFemtoEventReaderStandard* fstd = dynamic_cast<AliFemtoEventReaderStandard *> (fReader);
 387       if (fstd)
 388         {
 389           // Process the event with Kine information
 390           fstd->SetESDSource(fESD);
 391           if (mctruth) {
 392             fstd->SetStackSource(fStack);
 393             fstd->SetGenEventHeader(hdh);
 394             fstd->SetInputType(AliFemtoEventReaderStandard::kESDKine);
 395           }
 396           else
 397             fstd->SetInputType(AliFemtoEventReaderStandard::kESD);
 398           fManager->ProcessEvent();
 399         }
 400     }
 401
 402     // Post the output histogram list
 403     PostData(0, fOutputList);
 404   }
 405
 406   if (fAnalysisType==2) {
 407     if (!fAOD) {
 408       AliWarning("fAOD not available");
 409       return;
 410     }
 411
 412     // Get AOD
 413 //     AliAODInputHandler *aodH = dynamic_cast<AliAODInputHandler*> (AliAnalysisManager::GetAnalysisManager()->GetInputEventHandler());
 414
 415 //     if (!aodH) {
 416 //       AliWarning("Could not get AODInputHandler");
 417 //       return;
 418 //     }
 419 //     else {
 420
 421 //       fAOD = aodH->GetEvent();
 422 //     }
 423
 424     AliInfo(Form("Tracks in AOD: %d \n",fAOD->GetNumberOfTracks()));
 425
 426     if (fAOD->GetNumberOfTracks() > 0) {
 427       if (!fReader) {
 428         AliWarning("No AOD reader for AOD analysis! \n");
 429       }
 430       else {
 431         AliFemtoEventReaderAODChain* faodc = dynamic_cast<AliFemtoEventReaderAODChain *> (fReader);
 432
 433         if (faodc) {
 434           // Process the event
 435           faodc->SetAODSource(fAOD);
 436           fManager->ProcessEvent();
 437         }
 438         AliFemtoEventReaderStandard* fstd = dynamic_cast<AliFemtoEventReaderStandard *> (fReader);
 439
 440         if (fstd) {
 441           // Process the event
 442           fstd->SetAODSource(fAOD);
 443           fstd->SetInputType(AliFemtoEventReaderStandard::kAOD);
 444           fManager->ProcessEvent();
 445         }
 446       }
 447     }
 448
 449     // Post the output histogram list
 450     PostData(0, fOutputList);
 451   }
 452 }
 453
 454 //________________________________________________________________________
 455 void AliAnalysisTaskFemto::Terminate(Option_t *) {
 456   // Do the final processing
 457   if (fManager) {
 458     fManager->Finish();
 459   }
 460 }
 461 //________________________________________________________________________
 462 void AliAnalysisTaskFemto:: FinishTaskOutput() {
 463   // Do the final processing
 464   if (fManager) {
 465     fManager->Finish();
 466   }
 467 }
 468 //________________________________________________________________________
 469 void AliAnalysisTaskFemto::SetFemtoReaderESD(AliFemtoEventReaderESDChain *aReader)
 470 {
 471   AliInfo("Selecting Femto reader for ESD\n");
 472   fReader = aReader;
 473 }
 474 //________________________________________________________________________
 475 void AliAnalysisTaskFemto::SetFemtoReaderESDKine(AliFemtoEventReaderESDChainKine *aReader)
 476 {
 477   AliInfo("Selecting Femto reader for ESD with Kinematics information\n");
 478   fReader = aReader;
 479 }
 480 //________________________________________________________________________
 481 void AliAnalysisTaskFemto::SetFemtoReaderAOD(AliFemtoEventReaderAODChain *aReader)
 482 {
 483   AliInfo("Selecting Femto reader for AOD\n");
 484   fReader = aReader;
 485 }
 486 void AliAnalysisTaskFemto::SetFemtoReaderStandard(AliFemtoEventReaderStandard *aReader)
 487 {
 488   AliInfo("Selecting Standard all-purpose Femto reader\n");
 489   fReader = aReader;
 490 }
 491 void AliAnalysisTaskFemto::SetFemtoReaderKinematics(AliFemtoEventReaderKinematicsChain *aReader)
 492 {
 493   printf("Selecting Femto reader for Kinematics (Monte Carlo) information\n");
 494   fReader = aReader;
 495 }
 496 //________________________________________________________________________
 497 void AliAnalysisTaskFemto::SetFemtoManager(AliFemtoManager *aManager)
 498 {
 499   fManager = aManager;
 500   AliInfo(Form("Got reader %p\n", (void *) aManager->EventReader()));
 501   AliFemtoEventReaderESDChain     *tReaderESDChain     = dynamic_cast<AliFemtoEventReaderESDChain *> (aManager->EventReader());
 502   AliFemtoEventReaderESDChainKine *tReaderESDChainKine = dynamic_cast<AliFemtoEventReaderESDChainKine *> (aManager->EventReader());
 503   AliFemtoEventReaderAODChain     *tReaderAODChain     = dynamic_cast<AliFemtoEventReaderAODChain *> (aManager->EventReader());
 504   AliFemtoEventReaderStandard     *tReaderStandard     = dynamic_cast<AliFemtoEventReaderStandard *> (aManager->EventReader());
 505   AliFemtoEventReaderKinematicsChain *tReaderKineChain = dynamic_cast<AliFemtoEventReaderKinematicsChain *> (aManager->EventReader());
 506
 507  if ((!tReaderESDChain) && (!tReaderESDChainKine) && (!tReaderAODChain) && (!tReaderStandard) && (!tReaderKineChain)) {
 508     AliWarning("No AliFemto event reader created. Will not run femto analysis.\n");
 509     return;
 510   }
 511   if (tReaderESDChain) SetFemtoReaderESD(tReaderESDChain);
 512   if (tReaderESDChainKine) SetFemtoReaderESDKine(tReaderESDChainKine);
 513   if (tReaderAODChain) SetFemtoReaderAOD(tReaderAODChain);
 514   if (tReaderStandard) SetFemtoReaderStandard(tReaderStandard);
 515   if (tReaderKineChain) SetFemtoReaderKinematics(tReaderKineChain);
 516 }
 517