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