- //
- // Connect the Trees Kinematics and Hits for event # event
- // Set branch addresses
- //
-
- // Reset existing structures
- ResetHits();
- ResetTrackReferences();
- ResetDigits();
- ResetSDigits();
-
- // Delete Trees already connected
- if (fTreeH) { delete fTreeH; fTreeH = 0;}
- if (fTreeTR) { delete fTreeTR; fTreeTR = 0;}
- if (fTreeD) { delete fTreeD; fTreeD = 0;}
- if (fTreeR) { delete fTreeR; fTreeR = 0;}
- if (fTreeS) { delete fTreeS; fTreeS = 0;}
-
- // Create the particle stack
- if (fHeader) delete fHeader;
- fHeader = 0;
-
- // Get header from file
- if(fTreeE) {
- fTreeE->SetBranchAddress("Header", &fHeader);
-
- if (!fTreeE->GetEntry(event)) {
- Error("GetEvent","Cannot find event:%d\n",event);
- return -1;
- }
- }
- else {
- Error("GetEvent","Cannot find Header Tree (TE)\n");
- return -1;
- }
-
- // Get the stack from the header, set fStack to 0 if it
- // fails to get event
- //
- TFile *file = fTreeE->GetCurrentFile();
- char treeName[20];
-
- file->cd();
-
- if (fStack) delete fStack;
- fStack = fHeader->Stack();
- if (fStack) {
- if (!fStack->GetEvent(event)) fStack = 0;
- }
-
- // Get Hits Tree header from file
- sprintf(treeName,"TreeH%d",event);
- fTreeH = (TTree*)gDirectory->Get(treeName);
- if (!fTreeH) {
- //Error("GetEvent","cannot find Hits Tree for event:%d\n",event);
- Warning("GetEvent","cannot find Hits Tree for event:%d\n",event);
- }
-
- // Get TracReferences Tree header from file
- sprintf(treeName,"TreeTR%d",event);
- fTreeTR = (TTree*)gDirectory->Get(treeName);
- if (!fTreeTR) {
- //Error("GetEvent","cannot find TrackReferences Tree for event:%d\n",event);
- Warning("GetEvent","cannot find TrackReferences Tree for event:%d\n",event);
- }
-
- // get current file name and compare with names containing trees S,D,R
- TString curfilname=(TString)fTreeE->GetCurrentFile()->GetName();
- if(fTreeDFileName==curfilname)fTreeDFileName="";
- if(fTreeSFileName==curfilname)fTreeSFileName="";
- if(fTreeRFileName==curfilname)fTreeRFileName="";
- // Get Digits Tree header from file
- sprintf(treeName,"TreeD%d",event);
-
- if (!fTreeDFile && fTreeDFileName != "") {
- InitTreeFile("D",fTreeDFileName);
- }
- if (fTreeDFile) {
- fTreeD = (TTree*)fTreeDFile->Get(treeName);
- } else {
- fTreeD = (TTree*)file->Get(treeName);
- }
- if (!fTreeD) {
- // Warning("GetEvent","cannot find Digits Tree for event:%d\n",event);
- }
- if(fTreeDFileName != ""){
- if(fTreeDFileName==fTreeSFileName) {
- fTreeSFileName = "";
- fTreeSFile = fTreeDFile;
- }
- if(fTreeDFileName==fTreeRFileName) {
- fTreeRFileName = "";
- fTreeRFile = fTreeDFile;
- }
- }
-
- file->cd();
-
- // Get SDigits Tree header from file
- sprintf(treeName,"TreeS%d",event);
- if (!fTreeSFile && fTreeSFileName != "") {
- InitTreeFile("S",fTreeSFileName);
- }
- if (fTreeSFile) {
- fTreeS = (TTree*)fTreeSFile->Get(treeName);
- } else {
- fTreeS = (TTree*)gDirectory->Get(treeName);
- }
- if (!fTreeS) {
- // Warning("GetEvent","cannot find SDigits Tree for event:%d\n",event);
- }
-
- if(fTreeSFileName != ""){
- if(fTreeSFileName==fTreeRFileName){
- fTreeRFileName = "";
- fTreeRFile = fTreeSFile;
- }
- }
-
- file->cd();
-
- // Get Reconstruct Tree header from file
- sprintf(treeName,"TreeR%d",event);
- if (!fTreeRFile && fTreeRFileName != "") {
- InitTreeFile("R",fTreeRFileName);
- }
- if(fTreeRFile) {
- fTreeR = (TTree*)fTreeRFile->Get(treeName);
- } else {
- fTreeR = (TTree*)gDirectory->Get(treeName);
- }
- if (!fTreeR) {
- // printf("WARNING: cannot find Reconstructed Tree for event:%d\n",event);
- }
-
- file->cd();
-
- // Set Trees branch addresses
- TIter next(fModules);
- AliModule *detector;
- while((detector = (AliModule*)next())) {
- detector->SetTreeAddress();
- }
-
- fEvent=event; //MI change
-
- return fHeader->GetNtrack();
-}
-
-//_____________________________________________________________________________
-TGeometry *AliRun::GetGeometry()
-{
- //
- // Import Alice geometry from current file
- // Return pointer to geometry object
- //
- if (!fGeometry) fGeometry = (TGeometry*)gDirectory->Get("AliceGeom");
- //
- // Unlink and relink nodes in detectors
- // This is bad and there must be a better way...
- //
-
- TIter next(fModules);
- AliModule *detector;
- while((detector = (AliModule*)next())) {
- TList *dnodes=detector->Nodes();
- Int_t j;
- TNode *node, *node1;
- for ( j=0; j<dnodes->GetSize(); j++) {
- node = (TNode*) dnodes->At(j);
- node1 = fGeometry->GetNode(node->GetName());
- dnodes->Remove(node);
- dnodes->AddAt(node1,j);
- }
- }
- return fGeometry;
-}
-
-//_____________________________________________________________________________
-void AliRun::GetNextTrack(Int_t &mtrack, Int_t &ipart, Float_t *pmom,
- Float_t &e, Float_t *vpos, Float_t *polar,
- Float_t &tof)
-{
- // Delegate to stack
- //
- fStack->GetNextTrack(mtrack, ipart, pmom, e, vpos, polar, tof);
-}
-
-//_____________________________________________________________________________
-Int_t AliRun::GetPrimary(Int_t track) const
-{
- //
- // return number of primary that has generated track
- //
- return fStack->GetPrimary(track);
-}
-
-//_____________________________________________________________________________
-void AliRun::InitMC(const char *setup)
-{
- //
- // Initialize the Alice setup
- //
-
- if(fInitDone) {
- Warning("Init","Cannot initialise AliRun twice!\n");
- return;
- }
-
- gROOT->LoadMacro(setup);
- gInterpreter->ProcessLine(fConfigFunction.Data());
-
-
- gMC->DefineParticles(); //Create standard MC particles
-
- TObject *objfirst, *objlast;
-
- fNdets = fModules->GetLast()+1;
-
- //
- //=================Create Materials and geometry
- gMC->Init();
-
- // Added also after in case of interactive initialisation of modules
- fNdets = fModules->GetLast()+1;
-
- TIter next(fModules);
- AliModule *detector;
- while((detector = (AliModule*)next())) {
- detector->SetTreeAddress();
- objlast = gDirectory->GetList()->Last();
-
- // Add Detector histograms in Detector list of histograms
- if (objlast) objfirst = gDirectory->GetList()->After(objlast);
- else objfirst = gDirectory->GetList()->First();
- while (objfirst) {
- detector->Histograms()->Add(objfirst);
- objfirst = gDirectory->GetList()->After(objfirst);
- }
- }
- ReadTransPar(); //Read the cuts for all materials
-
- MediaTable(); //Build the special IMEDIA table
-
- //Initialise geometry deposition table
- fEventEnergy.Set(gMC->NofVolumes()+1);
- fSummEnergy.Set(gMC->NofVolumes()+1);
- fSum2Energy.Set(gMC->NofVolumes()+1);
-
- //Compute cross-sections
- gMC->BuildPhysics();
-
- //Write Geometry object to current file.
- fGeometry->Write();
-
- fInitDone = kTRUE;
-
- fMCQA = new AliMCQA(fNdets);
-
-// JCH note: the following line is useless, AliConfig instance is already
-// created in AliMC ctor. But it does not hurt.
- AliConfig::Instance();
- //
- // Save stuff at the beginning of the file to avoid file corruption
- Write();
-}
-
-//_____________________________________________________________________________
-void AliRun::MediaTable()
-{
- //
- // Built media table to get from the media number to
- // the detector id
- //
- Int_t kz, nz, idt, lz, i, k, ind;
- // Int_t ibeg;
- TObjArray &dets = *gAlice->Detectors();
- AliModule *det;
- //
- // For all detectors
- for (kz=0;kz<fNdets;kz++) {
- // If detector is defined
- if((det=(AliModule*) dets[kz])) {
- TArrayI &idtmed = *(det->GetIdtmed());
- for(nz=0;nz<100;nz++) {
- // Find max and min material number
- if((idt=idtmed[nz])) {
- det->LoMedium() = det->LoMedium() < idt ? det->LoMedium() : idt;
- det->HiMedium() = det->HiMedium() > idt ? det->HiMedium() : idt;
- }
- }
- if(det->LoMedium() > det->HiMedium()) {
- det->LoMedium() = 0;
- det->HiMedium() = 0;
- } else {
- if(det->HiMedium() > fImedia->GetSize()) {
- Error("MediaTable","Increase fImedia from %d to %d",
- fImedia->GetSize(),det->HiMedium());
- return;
- }
- // Tag all materials in rage as belonging to detector kz
- for(lz=det->LoMedium(); lz<= det->HiMedium(); lz++) {
- (*fImedia)[lz]=kz;
- }
- }
- }
- }
- //
- // Print summary table
- printf(" Traking media ranges:\n");
- for(i=0;i<(fNdets-1)/6+1;i++) {
- for(k=0;k< (6<fNdets-i*6?6:fNdets-i*6);k++) {
- ind=i*6+k;
- det=(AliModule*)dets[ind];
- if(det)
- printf(" %6s: %3d -> %3d;",det->GetName(),det->LoMedium(),
- det->HiMedium());
- else
- printf(" %6s: %3d -> %3d;","NULL",0,0);
- }
- printf("\n");
- }
-}
-
-//____________________________________________________________________________
-void AliRun::SetGenerator(AliGenerator *generator)
-{
- //
- // Load the event generator
- //
- if(!fGenerator) fGenerator = generator;
-}
-
-//____________________________________________________________________________
-void AliRun::ResetGenerator(AliGenerator *generator)
-{
- //
- // Load the event generator
- //
- if(fGenerator)
- if(generator)
- Warning("ResetGenerator","Replacing generator %s with %s\n",
- fGenerator->GetName(),generator->GetName());
- else
- Warning("ResetGenerator","Replacing generator %s with NULL\n",
- fGenerator->GetName());
- fGenerator = generator;
-}
-
-//____________________________________________________________________________
-void AliRun::SetTransPar(char *filename)
-{
- fTransParName = filename;
-}
-
-//____________________________________________________________________________
-void AliRun::SetBaseFile(char *filename)
-{
- fBaseFileName = filename;
-}
-
-//____________________________________________________________________________
-void AliRun::ReadTransPar()
-{
- //
- // Read filename to set the transport parameters
- //
-
-
- const Int_t kncuts=10;
- const Int_t knflags=11;
- const Int_t knpars=kncuts+knflags;
- const char kpars[knpars][7] = {"CUTGAM" ,"CUTELE","CUTNEU","CUTHAD","CUTMUO",
- "BCUTE","BCUTM","DCUTE","DCUTM","PPCUTM","ANNI",
- "BREM","COMP","DCAY","DRAY","HADR","LOSS",
- "MULS","PAIR","PHOT","RAYL"};
- char line[256];
- char detName[7];
- char* filtmp;
- Float_t cut[kncuts];
- Int_t flag[knflags];
- Int_t i, itmed, iret, ktmed, kz;
- FILE *lun;
- //
- // See whether the file is there
- filtmp=gSystem->ExpandPathName(fTransParName.Data());
- lun=fopen(filtmp,"r");
- delete [] filtmp;
- if(!lun) {
- Warning("ReadTransPar","File %s does not exist!\n",fTransParName.Data());
- return;
- }
- //
- if(fDebug) {
- printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
- printf(" *%59s\n","*");
- printf(" * Please check carefully what you are doing!%10s\n","*");
- printf(" *%59s\n","*");
- }
- //
- while(1) {
- // Initialise cuts and flags
- for(i=0;i<kncuts;i++) cut[i]=-99;
- for(i=0;i<knflags;i++) flag[i]=-99;
- itmed=0;
- for(i=0;i<256;i++) line[i]='\0';
- // Read up to the end of line excluded
- iret=fscanf(lun,"%[^\n]",line);
- if(iret<0) {
- //End of file
- fclose(lun);
- if(fDebug){
- printf(" *%59s\n","*");
- printf(" "); for(i=0;i<60;i++) printf("*"); printf("\n");
- }
- return;
- }
- // Read the end of line
- fscanf(lun,"%*c");
- if(!iret) continue;
- if(line[0]=='*') continue;
- // Read the numbers
- iret=sscanf(line,"%s %d %f %f %f %f %f %f %f %f %f %f %d %d %d %d %d %d %d %d %d %d %d",
- detName,&itmed,&cut[0],&cut[1],&cut[2],&cut[3],&cut[4],&cut[5],&cut[6],&cut[7],&cut[8],
- &cut[9],&flag[0],&flag[1],&flag[2],&flag[3],&flag[4],&flag[5],&flag[6],&flag[7],
- &flag[8],&flag[9],&flag[10]);
- if(!iret) continue;
- if(iret<0) {
- //reading error
- Warning("ReadTransPar","Error reading file %s\n",fTransParName.Data());
- continue;
- }
- // Check that the module exist
- AliModule *mod = GetModule(detName);
- if(mod) {
- // Get the array of media numbers
- TArrayI &idtmed = *mod->GetIdtmed();
- // Check that the tracking medium code is valid
- if(0<=itmed && itmed < 100) {
- ktmed=idtmed[itmed];
- if(!ktmed) {
- Warning("ReadTransPar","Invalid tracking medium code %d for %s\n",itmed,mod->GetName());
- continue;
- }
- // Set energy thresholds
- for(kz=0;kz<kncuts;kz++) {
- if(cut[kz]>=0) {
- if(fDebug) printf(" * %-6s set to %10.3E for tracking medium code %4d for %s\n",
- kpars[kz],cut[kz],itmed,mod->GetName());
- gMC->Gstpar(ktmed,kpars[kz],cut[kz]);
- }
- }
- // Set transport mechanisms
- for(kz=0;kz<knflags;kz++) {
- if(flag[kz]>=0) {
- if(fDebug) printf(" * %-6s set to %10d for tracking medium code %4d for %s\n",
- kpars[kncuts+kz],flag[kz],itmed,mod->GetName());
- gMC->Gstpar(ktmed,kpars[kncuts+kz],Float_t(flag[kz]));
- }
- }
- } else {
- Warning("ReadTransPar","Invalid medium code %d *\n",itmed);
- continue;
- }
- } else {
- if(fDebug) printf("%s::ReadTransParModule: %s not present\n",ClassName(),detName);
- continue;
- }
- }
-}
-
-
-//_____________________________________________________________________________
-void AliRun::MakeTree(Option_t *option, const char *file)
-{
- //
- // Create the ROOT trees
- // Loop on all detectors to create the Root branch (if any)
- //
-
- char hname[30];
- //
- // Analyse options
- const char *oK = strstr(option,"K");
- const char *oH = strstr(option,"H");
- const char *oTR = strstr(option,"T");
- const char *oE = strstr(option,"E");
- const char *oD = strstr(option,"D");
- const char *oR = strstr(option,"R");
- const char *oS = strstr(option,"S");
- //
-
- TDirectory *cwd = gDirectory;
-
- TBranch *branch = 0;
-
- if (oK) fStack->MakeTree(fEvent, file);
-
- if (oE && !fTreeE) {
- fTreeE = new TTree("TE","Header");
- // branch = fTreeE->Branch("Header", "AliHeader", &fHeader, 4000, 0);
- branch = fTreeE->Branch("Header", "AliHeader", &fHeader, 4000, 0);
- branch->SetAutoDelete(kFALSE);
- TFolder *folder = (TFolder *)gROOT->FindObjectAny("/Folders/RunMC/Event/Header");
- if (folder) folder->Add(fHeader);
-// branch = fTreeE->Branch("Stack","AliStack", &fStack, 4000, 0);
-// branch->SetAutoDelete(kFALSE);
-// if (folder) folder->Add(fStack);
- fTreeE->Write(0,TObject::kOverwrite);
- }
-
- if (file && branch) {
- char * outFile = new char[strlen(gAlice->GetBaseFile())+strlen(file)+2];
- sprintf(outFile,"%s/%s",GetBaseFile(),file);
- branch->SetFile(outFile);
- TIter next( branch->GetListOfBranches());
- while ((branch=(TBranch*)next())) {
- branch->SetFile(outFile);
- }
- if (GetDebug()>1)
- printf("* MakeBranch * Diverting Branch %s to file %s\n", branch->GetName(),file);
- cwd->cd();
- delete outFile;
- }
-
- if (oH && !fTreeH) {
- sprintf(hname,"TreeH%d",fEvent);
- fTreeH = new TTree(hname,"Hits");
- fTreeH->SetAutoSave(1000000000); //no autosave
- fTreeH->Write(0,TObject::kOverwrite);
- }